F/e  s/3 


AD-A118  496 
UNCLASSIFIED 


DELAWARE  UNZV  NEWARK  DEPT  OF  CXVlL  ENSINEERINS 

have  spectral  transformation  in  shallow  hater. 

JUN  62  Y  CHEN*  H  HANS  N00014-S1*K-0297 

CE-S2-23  NL 


DUC  FILE  COPY  AO  All  84 96 


% 


c 


WAVE  SPECTRAL  TRANSFORMATION 
IN  SHALLOW  WATER 

by 

Yun-Hal  Chen 
and 

Htlang  Wang 

t 


Technical  Report  No.  10 
Contract  No.  N00014-81-K-0297 

with  the  OFFICE  OF  NAVAL  RESEARCH  GEOGRAPHY  PROGRAMS 


DTIC 


Research  Report  CE-82-23 


OCEAN  ENGINEERING  PROGRAM 


PISTHIBimON  STATEMENT  A 

Approved  lot  public  release; 
Distiibution  Unlimited 


DEPARTMENT  OF  CIVIL  ENGINEERING 
UNIVERSITY  OF  DELAWARE 
NEWARK,  DELAWARE 
19711 


82  08  20  027 


wave  spectral  TRANSF0RMAT1<» 
XM  SHALLOW  WATER 


by 

Yun'Hal  Chen 
and 

Hsiang  Wang 


Accession  For 

NTIS  GRAJtl 
DTIC  TAB 
Unannounced 
Justification- 


By - - 

Dlstr 

Aval 

\  Dlat 

IL 

ibutlo 

labill 

lAvail 

Spec 

ty  Codes 

and/or 

ial 

June,  1982 

Technical  Report  Ho.  10 
Contract  Ho.  N00014-81-K-0297 
with  the  OPPICE  OF  NAVAL  RESEARCH  GEOGRAPHY  PROGRAMS 

Research  Report  CE”82”23 


OCEAN  ENGINEERING  PROGRAM 


DBPAR3WENT  OP  CIVIL  ENGINEERING 
UNIVERSITY  OP  DELAWARE 
NEWARK,  DELAWARE 
19711 


DISTRIBUTION  STATEMENT  A 

Approved  for  public  release) 
Distribution  Unlimited 


DO 


TABLE  OP  CONTENTS 


ACKNOWLEIXaiENTS . ill 

LIST  OP  FIGURES . vii 

KEY  TO  SY!4BOLS  AND  ABBREVIATIONS  .  X 

ABSTRACT . xiii 

CHAPTER .  I 

1.  INTRODUCTION  .  1 

X.l  Literaturs  Reviaw  .  2 

1.2  Objectives .  7 


2.  THEORY  OP  WAVE  SPECTRAL  TRANSFORMATION  .  .  10 

2.1  Description  of  a  Wave  Field  ....  10 

2.2  Conservation  of  Energy  Flux  in 

Wave  Spectral  Transformation  .  .  .  .  16 

2.3  Basic  Equations  for  Wave  Spectral 


Transformation . 17 

2.3.1  Wave  Kinematics . 17 

2.3.2  Wave  Dynamics . 22 

2.4  Specification  of  the  G  Functions  .  .  24 

2.4.1  Wave  Generation . 24 

2.4. 1.1  Phillips'  Resonant  Mechanism  .  25 

2.4. 1.2  Miles'  Instability  Mechanism  .  29 

2.4.2  Wave  Dissipation . 35 


V 


2. 4. 2.1  Bottom  Friction  .  35 

2. 4. 2. 2  Wav«  Breaking . 40 

3.  NUMERICAL  SCHEME  OF  STATIONARY  WAVE 

SPECTRAL  TRANSFORMATION  .  43 

3.1  Basic  Equations  for  Stationary  wave 

Spectral  Transformation  .  43 

3.2  Numerical  Procedures  .  45 

3.2.1  Wave  Kinematics  and  Boundary 

Conditions . 47 

3.2.2  Wave  Dynamics  and  Boundary 

Condition . 53 

3.3  Examples . 55 

3.3.1  Energy  Flux  Conservation  Case  .  .  57 

3.3.2  Bottom  Friction  Case . 60 

3.3.3  Local  Wind  Generation  Case  ...  63 

4.  NUMERICAL  SCHEME  OF  NON>STATIONARY  WAVE 

SPECTRAL  TRANSFORMATION  .  68 

4.1  Basic  Equations  for  Non-stationary 

Wave  Spectral  Transformation  ....  68 

4.2  Method  of  Solution . 70 

4.2.1  Two-Step  Lax-Wendroff  Scheme  .  .  70 

4.2.2  Stability  Analysis  .  73 

4. 2. 2.1  Homogeneous  Equation  ....  73 

4. 2. 2. 2  Non-Homogeneous  Equation  .  .  76 

4.2.3  Numerical  Dissipation  and 

Dispersion  Analysis  .  80 

4.3  Numerical  Analysis  for  Wave 

Conservation  Equation  .  86 


Vi 


4.3.1  Numerical  Procaduras  .  86 

4.3.2  Stablllt/  Analysis  .  90 

4.3.3  initial  Condition  and  Boundary 

Conditions . 91 

4.4  Nuraarisal  Analysis  for  Energy 

Transport  Equation  .  92 

4.4.1  Numerical  Procedures  .  92 

4.4.2  Stability  Analysis  .  94 

4.4.3  Initial  Condition  and  Boundary 

Conditions . 96 

4.5  Application  in  Monochromatic  Wave  .  .  96 

4.5.1  One-Dimensional  Casa . 98 

4.5.2  Two-Dimensional  Case . 105 

5.  RANDOM  WAVE  MODEL  FOR  NON-STATIONARY  WAVE 

SPECTRAL  TRANSFORMATION  .  112 

5.1  Swell  Wave  Model  . . 112 

5.2  Wind  Wave  Model . 116 

5.3  Hybrid  Wave  Model . 119 

5.4  Comparison  with  the  North  Sea 

Field  Data . 122 

6.  SUMMARY  AND  CONCLUSIONS . 135 

REFERRENCES . 142 


LIST  OF  FIGURES 


Figure 

1.1 

Schematic  View  of  Ocean  Environment 

Figxire 

2.1 

Schematic  View  of  Nearshore  Beach  Terminology 

Figure 

2.2 

Onshore  Wind  Blows  over  a  Sloping  Bottom 

Figure 

2.3 

Offshore  Wind  Blovrs  over  a  Sloping  Bottom 

Figure 

2.4 

Wind  Blows  over  a  Deep  Uniform  Bottom 

Figure 

3.1 

Full  Grid  and  Local  Grid  Description 

Figure 

3.2 

Schematic  View  of  Ocean  Environment 

Figure 

3.3 

Variation  of  Energy  Spectrum  in  Energy  Flux 
Conservation  Case  at  Boundary,  Location  5II  and 
Location  21X 

Figure 

3.4 

Variation  of  Energy  Spectrum  in  Energy  Flux 
Conservation  Case  at  Boundary,  Location  511  and 
Location  2II 

Figure 

3.5 

Percentage  Error  of  Spectral  Density  vs. 
Deviations  of  Wave  Directions 

Figure 

3.6 

Comparison  between  Energy  Spectrum  with  and 
without  Bottom  Friction  at  Location  5 II 

Figure 

3.7 

Comparison  between  Energy  Spectrum  with  and 
without  Bottom  Friction  at  Location  2II 

Figure 

3.8 

Energy  Spectrxjtm  of  Local  Wind  at  Location  SIX 
and  Location  2II 

Figure 

3.9 

Comparison  of  Energy  Spectrum  of  Local 

Wind  at  Location  5II  and  2II  with 
Pierson-Moskowitz  Spectrum 

Figure 

4.1 

Initial  Lax-Friedrichs  and  the  First  Leap-Frog 
Scheme  Solution  Steps 

j 


viii 


Fi9ur«  4.2 
Figur*  4.3 

Figure  4.4 

Figure  4.5 

Figtire  4.6 

Figure  4.7 

Figure  4.8 
Figure  4.9 

Figure  4.10 

Figure  4.11 

Figure  4.12 

Figure  4.13 

Figtore  4.14 

Figure  4.15 


General  Leap-Frog  Solution  Scheme 

Conqparison  between  Exact  and  Approximate 
Solutione  with  Hon-Homogeneous  Term  in  nth 
Time  Level  in  Leap-Frog  Scheme 

Ccrapariaon  between  Exact  and  Approximate 
Solutione  with  Hon-Bomogeneous  Term  in  (n-l)th 
Time  Level  in  Leap-Frog  Scheme 

Ccraparieon  of  Exact  and  Numerical  Solution 
for  r«l/2 

Dissipation  and  Dispersion  for  the 
Lax-Friedrichs  Scheme 

Dissipation  and  Dispersion  for  the 
Leap-Frog  Scheme 

Schematic  View  of  a  Two  Dimensional  Case 

Variation  of  Wave  Energy  Decreasing  with 
T»10sec*  8«0‘*at  Boundary,  Location  5III  and 
Location  2111 

Variation  of  Wave  Energy  Decreasing  with 
T«2.5sec,  9 at  Boundary,  Location  5111  and 
Location  2111 

Variation  of  Wave  Energy  Increasing  with 
T«10sec,  SbO**  at  Boundary,  Location  5111  and 
Location  2111 

Variation  of  Wave  Energy  Decreasing  with 
T«10sec,  9»20°at  Boundary,  Location  5111  and 
Location  2111 

Variation  of  Wave  Angle  Decreasing  with 
TslOsec,  9*20*’ at  Location  5III  and  Location 
2111 

Variation  of  Wave  Angle  Decreasing  with 
T*10»ec,  9*-20*at  Location  5111  and  Location 
2111 

Variation  of  Wave  Energy  Increasing  with 
T^lOsec,  9*20" at  Boundary,  Location  5111  and 
Location  2111 


Figure  4.16  Variation  of  Wave  Angle  Increasing  with 

T»10aec<  O^icf  at  Location  5III  and  Location 
2III 

Figure  5.1  Schematic  View  of  Ocean  Environment 

Figure  5.2  Variation  of  Energy  Spectrum  Decreasing  due  to 
Swell  at  Boundary,  Location  5II  and  Location 
2II 

Figure  5.3  Variation  of  Energy  Spectrxjm  Increasing  due  to 
Swell  at  Boundary,  Location  5II  and  location 
2II 

Figure  5.4  Time  Changes  of  Wave  Energy  Spectrum  in  the 
Decrease  of  Swell  of  a  Typhoon  (After  Ijima, 
1957) 

Figure  5.5  Variation  of  Energy  Spectrum  Increasing  due  to 
Wind  Generation  at  Location  5 II  and  Location 
2II 

Figure  5.6  Variation  of  Energy  Spectrum  due  to  Decrease  of 
Swell  and  Increase  of  Wind  Wave  at  Location  5II 
and  Location  2II 

Figure  5.7  Frequency  Spectrum  Indicating  Wind  Sea  and 

Swell  Peaks  (After  Long  and  Hasselmeunn,  1979} 

Figure  5.8  MARSEN  Wave  Measurement  Stations 

Figure  5.9  Nearshore  Field  Instriunent  Arrangement 

Figure  5.10  Wind  Speed  and  Direction  on  September  24,  1979 

Figure  5.11  Energy  Spectrum  and  Mean  Wave  Direction  at  Deep 
Water  at  12:14,  September  24,  1979 

Figure  5.12  Grid  System  used  in  the  Computation 


Figure  5.13  Comparison  between  Swell  Wave  Model  from  13:14, 
Wind  Wave  Model  from  8:00  and  Field  Data  at 
14:00  at  Location  122Sm 

Figure  5.14  Comparison  between  Hybrid  Wave  Model  and  Field 
Data  at  14:00  at  Location  1225m 


KEY  TO  SYMBOLS  AND  ABBREVIATIONS 


A(f)  Fraquency  spectral  density  function 

A(f,9)  Directional  frequency  spectral  density  function 

a  Wave  amplitude 

B  Coefficient  of  Phillips*  equilibrium  spectrum 

B(r,a)  Cross  spectral  density  function 

C  Wave  celerity 

Cg  Wave  group  velocity 

Cf  .  Bottom  friction  coefficient 

Co  Co  spectrum 

5  ..Total  energy  per  unit  surface  area 

F(lc)  Wave  number  spectral  density  function 

F(k,9)  Directional  wavs  number  spectral  density  function 
f  frequency 

6  Wave  generation  or  dissipation  function 

g  Acceleration  of  gravity 

h  Mean  water  depth 

i  Mesh  point  in  x-axis 

j  Mesh  point  in  y-axis 

k  Wave  number 

Wave  length 


X 


L 


M 

N 

n 

P 

P(lc,tf) 

q 

R 

r 

Sxx 

Sxy 

Syy 

T 

t 

At 

u 

Uc 


V 

w 

X 

AX 

y 

Ay 


xi 

Humbar  of  space  steps  in  x^axis 
Number  of  apace  steps  in  y-axis 
Mesh  point  in  time  axis 
Pressure 

Three  dimensional  spectrum  of  the  random 
atomospheric  pressure  fluctuation 

Quadrature  spectrum 

Co-variance  function 

Vector  space 

Radiation  stress  in  xx  direction 
Radiation  stress  in  xy  direction 
Radiation  stress  in  yy  direction 
Wave  period 
Time 

Time  step 

Current  velocity  in  x  direction 

Convection  velocity  of  the  the  turbulent 
eddies  with  wave  number  k 

Horizontal  water  partical  velocity  at  the  bottom 
Current  velocity  in  y  direction 
Wind  velocity 

Horizontal  axis  toward  offshore  direction 
space  step  in  x-axis 

Horizontal  axis  toward  alongshore  direction 
space  step  in  y-axis 


Vertical  axis  toward  upward 

Wave  direction 

Density  of  air 

Density  of  water 

Wave  angular  frequency 

Wave  angular  frequency 

for  stationary  reference  frame 

Phase  function 

Angle  between  wind  and  wave 

Phillips '  resonant  spectrum 

Miles'  instability  spectrum 

Coefficient  of  Pier son -Moskowitz  spectrum 

Coefficient  of  Pierson-Moskowitz  spectrum 

Total  bottom  friction  of  the  wave  field 

Water  surface  displacement 

Time  lag 

Bottom  friction 

Angular  Distribution 

Coefficient  for  breaking 

Representation  shown  in  Eq.  (  2 >44) 

Kitaigorodskii ' s  modification  factor 


AaSTR\CT 


The  pucpose  of  this  study  is  to  develop  s 
numerical  method  to  compute  shallow  water  wave  spectral 
transformation  both  for  stationary  and  non-s tationary 
cases  . 


A  numerical  model  for  stationary  wave  spectral 
transformation  has  been  developed  by  Shiau  and  Wang 
(1977)  and  further  amplified  by  Wang  and  Yang  (1991). 
The  present  work  extends  the  model  to  add  bottom 
friction  and  local  wind  generation.  This  model  is 
applicable  to  wave  spectral  transformation  under 
stationary  meterological  conditions.  If  the 
meterological  condition  is  non-s tationary  in  the 
generating  area  or  the  local  wind  effect  in  the  domain 
of  consideration  can  no  longer  be  neglected,  such  as 
when  a  hurricane  or  storm  sweeps  through  offshore  of  a 
coastal  region,  the  wave  spectral  transformation  in 
shallow  water  should  also  be  treated  as  non-s tationary . 
A  non-s tationary  wave  spectral  transformation  model  is 
thus  developed  to  handle  these  cases  . 
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xiv 

The  non-atationary  model  consists  of  three 
sub-models;  dealing  respectively,  with  1)  swell 
transformation  due  to  offshore  wind  variation  without 
local  wind  generation,  2)  wind  wave  transformation  due 
to  local  wind  generation,  3}  hybrid  condition  that 
combines  the  swell  and  the  wind  wave  transformations. 
The  numerical  results  were  compared  with  North  Sea 
field  data. 
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CHAPTER  1 


IMTRODUCTIOM 


The  need  for  a  more  realistic  representation 
and  accurate  prediction  to  shallow  water  environment  is 
increasing  for  both  scientific  and  engineering 
purposes.  It  is  also  an  accepted  opinion  that  wave 
spectrum  utilizing  the  ensemble  of  random  surface 
oscillations  provide  a  more  realistic,  and  perhaps  more 
accurate  representation  of  ocean  waves  than  a 
monochromatic  wave  train  of  single  frequency.  In  the 
past,  considerable  efforts  have  been  devoted  to  the 
develoEKnent  of  deepwater  wave  spectra,  to  the  extent 
that  reasonable  confidence  can  be  placed  in  its 
applicability  to  describe  wind  generated  ocean  waves. 
To  establish  the  wave  climate  in  the  coastal  zone,  as 
is  often  necessary  for  scientific  and  engineering 
purposes,  the  deepwater  wave  environment  needs  to  be 
transformed  into  shallow  water.  It  is  in  this  area 
that  both  knowledge  and  effort  are  conspicuously 
lacking,  although,  odd  as  it  may  seem,  the 
transformation  of  monochromatic  wave  trains  is  a 
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vrall-studied  topic. 

l.I  Literatara  Review 

There  exist  a  few  attempts  in  the  development 
of  theoretical  concepts  and  numerical  techniques 
concerning  wave  spectral  transformation  from  deep  to 
shallow  water . 

Longuet-Higgins  (1956,1957)  considered  the 
refraction  of  a  wave  spectrum  having  just  two 
components.  In  one  case,  the  two  components  consisted 
of  two  regular  wave  trains  of  the  same  wave  number  and 
travelling  in  slightly  different  directions.  The 
effects  of  refraction  on  each  other  tend  to  increase 
the  crest  lengths,  and  to  decrease  the  wave  heights  and 
the  angle  between  the  two  components.  Another  case 
considered  by  Longuet-Higgins  was  a  two  component 
spectrum  where  the  components  were  of  slightly 
different  wave  numbers  and  both  travelling  in  the  same 
direction  in  deep  water.  In  this  case,  the  wave 
pattern  which  originally  was  long  crested  becomes  short 
crested  due  to  the  different  refraction  of  the  two 


components . 
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Bassd  on  the  radiative  transfer  equation. 


Barnett 

(1968) 

developed  a 

method 

to 

predict  the 

two-dimens ional 

wind 

wave 

spectrum 

in 

the  Nor  th 

Atlantic 

Ocean. 

The 

model 

takes 

account  of  wave 

generation  by  both  resonance  and  instability  mechanisms 
and  vave  dissipation  by  breaking.  However,  the  shallow 
water  effects  such  as  refraction,  bottom  dissipation 
and  the  non-s tationary  wave  conservation  phenomenon  are 
not  involved.  Therefore,  the  model  is  applicable  to 
deep  water  only. 

Based  on  the  conservation  of  energy  flux, 
Karlsson  (1969)  developed  a  method  to  compute  spectral 
transformation  for  refraction  over  parallel  bottom 
contours.  The  governing  aquation  for  the  distribution 
of  the  continuous  directional  wave  spectrum  under 
steady  state  conditions  in  water  of  any  depth  was 
derived  and  its  solution  by  finite  difference 
techniques  was  presented. 

Collins  (1972)  extended  the  work  to  the 
inclusion  of  bottom  frictional  effects  and  mentioned 
the  application  to  an  irregular  bottom  topography. 
However,  his  numerical  scheme  traced  wave  energy  along 
wave  rays,  which  made  the  computational  procedure  quite 


impractical  to  datarmina  shallow  watar  wava  spactra  at 
dasignatad  locations  in  an  araa  of  irragular  offshora 
bottcxn  configuration. 

Basad  on  a  geometrical-optical  approximation, 
Krasitskiy  (1974)  has  derived  explicit  analytical 
solutions  for  the  spectral  transformation  over 
two-dimems ional  parallel  bottom  contours.  The  solution 
he  achieved  took  account  of  refraction  and  was  applied 
to  problems  of  spectral  transformation  due  to 
diffractions . 

Hasselmann  (1976)  developed  a  parametric  model 
within  which  non-linear  interactions  are  implicitly 
taken  into  account  through  the  assumed  shape  of  the 
spectrum  represented  by  a  limited  number  of  parameters . 
Although  the  parametr ization  of  the  non-linear 
wave-wave  Interactions  had  been  considered,  shallow 
water  phenomenon  due  to  the  interaction  of  surface 
waves  with  the  bottom  was  not  taken  into  account. 

Shiau  and  Wang  (1977)  developed  a  numerical 
model  to  compute  the  stationary  wava  spectral 
transformation  over  irregular  bottom  topographies.  It 
is  based  on  the  assumption  that  wava  energy  associated 
with  a  narrow  frequency  band  stays  within  the  band  on 
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refraction  when  waves  propagate  from  deep  into  shallow 
water . 

Longuet-Higgins  (1957),  assuming  a  horizontal 
bottom,  obtained  an  approximate  solution  of 
monochromatic  wave  In  shoaling  water .  Yang  and  Chen 
(1979)  further  Improved  the  solution  by  incorporating 
wave  statistics.  They  showed  that  the  transformation 
formula  of  random  wave  spectrum,  based  upon  the 
asymptotic  solution  of  Friedrichs'  (1948)  small  beach 
slope  case,  was  consistent  with  the  approximation 
theory  by  Longuet-Hlggins •  In  addition,  they  extended 
Stoker's  (1947)  sloping  bottom  solution  from 
deterministic  wave  to  a  random  wave  field.  Their 
solution  was  based  on  the  consideration  of  stationary 
case  only;  the  non-s tatlonary  case  was  not  considered. 

Gunther  (1979)  developed  a  numerical  wave 
prediction  model  incorporating  a  parametrical  wind-sea 
model  and  a  characteristic  swell  model.  However,  the 
model  doesn't  include  the  shallow  water  effects  such  as 
refraction,  bottom  dissipation  and  also  can  not  handle 
the  changing  wind  effect. 
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Wang  and  Yang  (1981)  compared  field  results 
with  the  numerical  computations  for  energy 
transformation  in  shallow  water,  based  on  the  numerical 
model  developed  by  Shiau  and  Wang  (1977).  The  results 
suggested  that,  prior  to  wave  breaking,  the 
transformation  of  wave  components  is  mainly  influenced 
by  shoaling  and  refraction.  The  bottom  frictional 
effect  could  also  be  Important  in  the  energy-containing 
range.  After  breaking,  turbulent  generation  due  to 
wave  instability  seems  to  be  the  dominant  mechanism  of 
energy  dissipation. 

Based  on  the  ray  technique,  Cavaleri  and 
Riszoli  (1981)  developed  a  wind  wave  prediction  model 
in  shallow  water  which  included  wave  refraction, 
shoaling,  generation,  and  dissipation  (breaking  and 
bottom  friction) .  However,  without  considering  the 
wave  conservation  equation,  the  swell  condition  can  not 
be  handled. 

Vincent  (1982)  attempted  to  summarize  the  start 
state  of  art  on  shallow-water  wave  modeling  by 
examining  various  factors  that  have  been  and  that 
should  be  included  in  the  model.  The  review  was 
restricted  to  models  using  wave-ray  methods.  He 


observed  that  a  complata  shallow  watar  wava  modal 
should  Includa  1)  propagation  with  rafraction  and 
shoaling,  2)  sourca  machanism  dua  to  atmospharaic 
input,  3)  sink  machanism  including  bottom  dissipation, 
bottom  parcolation  and  wava  braaking  and  4)  non-linaar 
wava-wava  intar actions . 

1.2  Objactivas 

A  major  motivation  of  tha  prasant  study  is  to 
davalop  an  accurata  mathod  to  pradict  wava  spactral 
transformation  in  shallow  watar  anvironmant.  Tha 
ragion  of  intarast  is  dapictad  in  Figura  1.1  which 
bagins  from  tha  point  shoraward  of  tha  ganaratlng  araa 
to  tha  point  bafora  braaking.  Two  diffarant  approachas 
ara  us ad  in  tha  prasant  work.  Tha  first  mathod  is 
basad  on  Noda's  (1974)  ralaxation  finita  diffaranca 
schama  to  solva  tha  stationary  wava  spactral 
transformation  in  shallow  watar  with  tha  inclusion  of 
bottom  dissipation  and  local  wind  ganarabion  affacts . 
Tha  sacond  mathod  is  using  two-stap  Lax-Wandroff  schama 
to  solva  tha  non-s tat ionary  casa  of  wava  spactral 
transformation  in  shallow  watar. 


Figure  1.1  Schematic  View  of  Ocean  Environment 


The  stationary  wave  spectral  transformation 
method  is  applicable  to  the  energy  spectral 
transformation  under  stationary  meterological 
condition.  However,  if  the  meterological  condition  is 
non-s tat ionary  in  the  generating  area  and/or  in  the 
wave  transformation  area,  the  non-s tationary  wave 
spectral  transformation  method  should  be  applied. 

For  the  non-s tationary  wave  spectral 
transformation  method,  it  consists  of  three  sub-models 
:  L)  the  swell  model  which  handles  non-s tationary 
swell  transformation  without  local  wind  generation,  2) 
the  wind  wave  model  which  deals  with  the  wave 
transformation  due  to  wind  generation  effect,  3)  the 
hybrid  wave  model  that  combines  the  swell  wave  model 
and  the  wind  wave  model  to  handle  the  mixed  conditions . 
A  set  of  field  data  recently  collected  along  North  Sea 
was  compared  with  the  numerical  results. 


CHAPTER  2 


THEORY  OF  WAVE  SPECTRAL  TRANSFORMATION 

As  ocean  surface  waves  propagate  from  deep  to 
shallow  water,  significant  changes  occur  in  their 
characteristics.  These  changes  include  both  wave 
kinematic  and  dynsunic  properties.  Kinematically,  the 
wave  form  obeys  wave  conservation,  dispersion,  relation 
and  wave  refraction  phenomena.  Dynamically,  the  wave 
energy,  hence  the  wave  height,  changes  due  to  wave 
refraction,  wave  shoaling,  wave  generation,  and  wave 
dissipation.  Thus,  in  considering  wave  transformation 
in  shallow  water,  both  wave  kinematic  and  wave  dynamic 
properties  should  conform  to  governing  equations  which 
are  developed  in  this  chapter. 

2.1  Description  of  A  Wave  Field 

The  waves  appearing  on  the  surface  of  the  sea 
are  almost  always  random  in  the  sense  that  the  detailed 
configuration  of  the  surface  varies  in  an  irregular 
manner  in  both  space  and  time.  Only  the  various 
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statiatical  meaatarea  of  the  motion  can  be  regarded  aa 
significant  observationally  or  predictable 

theoretically.  Of  these,  the  fundamental  measure  is 
the  joint  probability  density  of  the  variables 
concerned.  In  a  homogeneous  and  stationary  wave  field, 
all  the  joint  probability  densities  are  spatially  and 
time  invariant  except  the  addition  of  a  constant  space 
vector  or  a  constant  time  scalar.  If  7i  ,  7/i  ere 

the  ocean  surface  variations  at  n  points,  then  the 

probability  density  function  of  7<  cen  be  expressed  as 

P(7..7a  ••d7a  -  1  i2J) 

where 

P ( 7. , *&.  •  .Ti)  is  the  probability  density  of  the 
surface  displacement. 

P(7i«  V  • 'TiJ  dud Tr  ••  (37m  represents  the  probability 
that  the  surface  displacements  at  the  points 

(S|  ,t|  )  •  .  .(XM,tM)  • 

i.a.  the  total  probability  is  equal  to  1. 

It  is  usually  convenient  to  measure  the 

displacement  'J  from  the  mean  free  surface  level,  so 
that  the  first  moment  (mean)  of  the  probability  density 
P( 7r)  becomes 

7r  “  7r  Pi7r)<!l7r  -  0 


r 

1 


i 
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for  all  r 

where  the  overbar  is  used  to  denote  the  ensemble 
average . 


The  second  moment  (covariance)  of  the 
probability  density  P(T..lt)  becomes 

R(5.r;t.T)  -  j  |  ^  7. )  flil.  ^7* 

*»  t  7  t  r  )  (2.3) 

where  the  points  1,2  are  taken  as  (x,t),  (x-<>r,t-fr  ), 

respectively.  In  particular,  the  mean  square  surface 
displacement  can  be  expressed  by  the  covariance  when 
both  and  X  *0 ,  i . e . , 


—  R  U.o  ;  t,®)  —  J  7/PC7.)<^7, 


U.4) 


The  wave  spectrvon  k,  »•  )  and  the  covariance 
R(r,T)  form  a  Fourier  transform  pair: 

OC  U.®-)  =»  Ur)  *  Jj  R(r.Oexp(-ilh.f-o-t)j 

(.2.5) 


and 


R  (r,  r)  «  U  0CU.®')exp[  i(k.r-et)]  dkdo. 
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Equations  (2.4)  and  (2.6)  lead  to 

1‘ =  a  (k,^)dkd»- =  E  (»7) 

where 

k  is  wave  number  vector,  |KI  *  k  * 

2E/T,  is  wave  angular  frequency( rad/s ec) , 

L  is  wave  length, 

T  is  wave  period, 

E  ■  7*  *  l/2a*  (except  for  a  factor  ^g)  is  the 

total  energy  per  unit  surface 
area, 

a  is  wave  amplitude. 

Therefore,  0(.(k,^)  can  be  interpretated  as  the  spectral 
density  function  per  unit  surface  area  of  wave  number 
and  angular  frequency. 

Reduced  spectral  density  functions  can  be 
obtained  from  OC(k,^)  by  integration  over  or  over  k  . 
The  wave  number  spectral  density  function  is  given  by 


FCk)  =  ' 

'■  o» 

i  a  (  K  .0^) 

(2.«) 

The  angular  frequency  spectral  density 

function  is 

given  by 

r 

A^o^)  =  J 

1  oc  (k.> )  d  K 

{2.9) 
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The  spectral  density  function  F(k)  and  A(^)  can 
be  considered  as  a  descriptor  of  the  sea  in  space  (x,y) 
and  time  t,  respectively.  F(k),for  instance,  is  a  five 
dimensional  variables  F(k.  ,ks ,x,y, t) .  The  angular 
frequency  spectral  density  function  is  also  a  five 
dimensional  variables,  A(0>,QjX,y, t)  with  Q  being  the 
directional  angle.  Alternate  forms  of  the  energy 
density  function  can  be  shown  as  follows: 

In  directional  wave-number  spectrum, 

F  (k  .e)  -  K  F 


where 


k  .  »  k  cos  8  ontk  ka  “  k  SJo  9 


In  directional  angular  frequency  spectrum 

A  (0^,9)  -  7^  F  (k.e) 


where  Cg 


dk 


is  the  group  velocity. 
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The  directional  frequency  spectrum  is 


A  (f.9)  -  an.  A  lo.  6) 


(2.12) 


where  f  is  wave  frequency  and0>«2]lf. 
The  one  dimensional  wave  spectrum  is: 

A  (f )  -  A  (f.e)  <ie 


(2.13) 


Therefore,  the  transformation  between  wave  number 
spectrum  and  directional  frequency  spectrum  becomes 


F(n..k.)=  F(k) --^*7  A  (f.e)  'i  'so) 


It  is  also  noted  that  by  definition. 


7*  *  f  f  F  (  k, ,  k,)  dk.  dkj 


f  f  A(f.9  )  df  de 


I2.l3k) 


Where  7**  is  the  value  of  the  mean  square  sea  surface 
fluctuation.  The  directional  frequency  spectrum  will 
be  used  throughout  this  text  since  it  is  easy  to 


16 


calculate.  • 


2.2  Conservation  of  Energy  Flux  in  Wave  Spectral 
Trans  formation 


Let  subscript  0  denote  the  reference  state 
(deepwater  condition  is  usually  conveniently  selected), 
the  conservation  aquation  of  energy  flux  transformation 
is  established.  Prom  equation  (2.13),  we  have 


A(f.e) 


fK  F(Kli<IKl 

(KF^k)dK), 


(A(f  9)df]. 


(2.14) 


For  steady  state  condition,  the  frequency  f  is 
invariant  in  space,  the  above  equation  simplifies  to 


ACfe) 


(k  FC(c)otKl 

(kF(R)<<KL 


A.Cf.e) 


since  F(k)d]c  approximates  the  incremental  of  wave 
energy  over  the  wave  number  band  dk,  it  is  plausible  to 
write 


dE  -  FUOdK 


C2.»t ) 


Substituting  aquation  (2.16)  into  equation  (2.15) 
results  in 


(KdE], 


A.(f  e; 


K3.I7) 
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This  transformat.ion  aquation  is  not  valid  whan 
wava  ganaration  and  dissipation  ara  involvad.  Also,  it 
cannot  ba  appliad  to  unstaady  stata  case  because  f  is 
no  longer  space  invariant.  More  details  will  ba  shown 
in  the  following  sections. 

2.3  Basic  Equations  for  Wave  Spectral  Transformation 

The  basic  equations  for  wave  spectral 
transformation  consist  two  parts:  one  is  wave 

kinematics  and  the  other  is  wave  dyneunics.  In  the 
derivations,  wave  reflection  and  non-linear  wave-wave 
interactions  are  neglected. 

2.3.1  Wave  .Kinematics 

Starting  with  a  progressive  linear  gravity 
wave,  the  free  surface  can  be  written  as 

'(U.y.t)**  a  (ir.y.t)  &»5  ^ ^ 

where  a  is  the  wave  amplitude  and  ^  «k>x-»t  is  phase 
function . 

A  wave  number  vector  can  be  defined  as 


k  - 


U.I8) 
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and  a  wav#  fraqucncy  3 c alar  can  ba  daflnad  as 


Assuming  ^(x,y,t)  Is  continuous,  then  the  order  of 
differentiation  are  interchangeable; 


Substituting  Equation  (2.13)  and  (2.19)  into  the  above 
expression,  it  is  found  that 


^  +  yo.  =  0 


U.2O 


which  is  the  classical  conservation  of  wave  equation. 
If  no  new  waves  are  being  created  by  a  local 
disturbance,  the  rate  of  change  of  wave  number  is 
balanced  by  the  convergence  of  the  frequency,  the 
number  of  wave  crests  passing  a  point  per  unit  time,  or 
the  flux  of  waves .  In  a  random  wave  field  of  linearly 
superposed  waves.  Equation  (2.21)  holds  for  each 
Fourier  component  but  the  simple  geometrical 
interpretation  is  not  possible.  If  the  wave  field  is 
constant  in  time,  thenT^“0,  or,  the  wave  period  does 
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not  Chang*  with  spac*.  It  ramaln*  constant  mvmn  as  th* 
wat*r  d*pth  changes. 

In  th*  presence  of  a  steady  current  u«ut*»>v5#  it 
can  be  shown  that  the  scalar  frequency  with  respect  to 
a  stationary  reference  frame  is 

OL  —  O  +  K-U 


The  wave  frequency  O'  with  respect  to  a  moving  reference 
frame  is  given  by  the  dispersion  relation. 


^  k  'taeh  kh 


If  it  is  also  assumed  that  the  wave  number  field 
changes  slowly  with  time  then  from  Equation (2. 21) 


Vl^+k  u)  **  0 


or 


^  +  k .  U  »■  Censtont 


(2.24) 


This  constant  can  be  evaluated  for  the  case  where  u>bO 
in  which  case  0’»2K/T  where  T  is  the  wave  period. 
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Equation  (2.24)  bacomes 


t  K  -  U 


or 


(SktonhkK)^  +  Ukcase  +  \/kJ.«e--^  C2.25) 


for  the  coordinate  system  shown  in  Figure  2.1.  Here  U 
and  V  are  current  velocity  components  in  the  x  and  y 
directions,  respectively;  is  the  wave  approach  angle 
measured  from  the  positive  x-axis  . 

Using  the  mathematical  property  that  the  curl 
of  a  gradient  is  identically  zero,  it  is  shown  that 

V  X  V  ^  ^  (2.  at> 

which  implies  that 

^XK«*0  (2.27) 


This  equation  states  that  the  wave  number  vector 
irrotational .  Using  the  coordinate  in  Figure 
Equation  (2.27)  can  be  further  decomposed  to 

^(KcosO?  __  <^(KS.n9)  ^  Q 

ay 


is 


2.1, 


\2.28) 
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For  a  shora  Ilna  whara  tha  alongshora  variations  in  tha 
y  diraction  of  all  variablas  ara  zaro,  that  is,  thara 
ara  straight  and  parallal  off shora  contours,  this 
aquation  yialda  Snail's  law. 

With  lcno%m  U,  V,  and  h,  tha  thraa  unknoim  wava 
nunbar  k,  angular  fraquancyo»,  and  wava  diraction  9  can 
ba  solvad  through  Equations  (2.21),  (2 >25)  and  (2.28). 
Howavar,  the  kinematic  conditions  can  not  represent 
energy  spectrum  unless  tha  dynamic  condition  is  also 
considered . 

2 *3. 2  Wava  Dynamics 

To  determine  the  variation  of  A,  the  assumption 
is  made  that  the  wave  energy  associated  with  a  certain 
frequency  band  stays  within  this  band  so  that  the 
linear  superposition  can  be  used.  This  assumption 
denies  energy  transfer  among  different  frequencies  and 
is  compatible  only  when  non-linear  effects  are 
negligible. 

Gelci,  et  al.  (1956),  Hasselmann  (1960),  and 
Groves  and  Melcer  (1961)  have  all  independently 
proposed  a  transport  equation  which  takes  into  account 
wave  propagation  as  well  as  the  processes  of  generation 
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and  dissipation.  Howavar,  tha  shoaling  of  wavas  is 
also  affastad  by  tha  intaractlons  of  wavas  and 
currants.  Tha  anargy  transport  aquation  prasantad  by 
Xonguat -Higgins  and  Sta%irart  (1960,1961)  for  tha 
fluctuating  motion  of  wavas  with  a  suparimposad  currant 
systam  is 


Ss)f  yy 


-  2 

i«l 

whaca 

A  ■  A(f,9,x,y,t)  is  directional  frequency  spectrum 
G  *  G(f,9/X,y,t)  is  a  general  representation  of 
all  processes  which  are  adding  or  subtracting 
energy  from  the  spectrum 

Cg  is  the  group  velocity  expressed  as 

_  1  ..  2kh  . 

C9  -  nC  -  T<^*JJS;ii53 

C"(  is  the  wave  celerity 

are  the  radiation  stresses  such  that: 

Sxx  «  \(  ( 2n-L/2  )cos*9  ( n-1/2  )3  in*9 )  ■  \9>xx 
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Syy  -  A(  (2n-l/2)«in*0  +  (n-l/2)coa*a  )  -  A<ryy 
Sxy  ■  Syx  •  Ancoa0ain0  •PJmy 
Tha  Equationa  (2.21) >  (2.25),  (2.28)  and  (2.29)  ara  tba 
baaic  aquatlona  for  wava  apact.ral  tr ana  formation . 

2.4  Spaciflcation  of  tha  G-Functlona 

In  Equation  (2.29),  G-Function  rapreaanta  all 
procaaaaa  which  ara  adding  energy  to  or  aubtracting 
energy  from  tha  apectrum,  i.e.,  aource  or  aink 
function.  Tha  proceaaea  can  be  categorized  into  linear 
and  nonlinear  mechanism.  The  detaila  are  as  follows. 

2.4.1  Wava  Generation 

Based  on  the  dynamics  of  the  surface  boundary 
layer,  both  the  surface  pressure  and  shear  stress 
variations  influence  wave  growth.  One  can  see 
intuitively  that  the  pressure  fluctuations  on  the 
surface  are  of  two  kinds:  those  produced  by  the 

turbulent  eddies  in  the  wind  and  those  induced  by  the 
air  flow  over  the  irregular  water  surface.  At  present, 
there  are  two  prominent  postulations  on  wind  wave 
generations . 
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2.4.1  Phillips'  Rssonant  Mschanism 

In  tha  rssonanca  thaory  proposad  by  Phillips 
(19S7),  tha  turbulant  atmospharic  pressura  fluctuations 
ara  saan  as  prims  wava  ganarators.  It  is  first  assumed 
that  tha  ocaan  and  tha  atmosphere  are  an  uncoupled 
system.  Provided  tha  turbulent  pressure  spectrum  is 
reasonably  isotropic,  the  main  condition  for  resonance 
is  that  Uccos  S  *  C(  O' ) ,  where  C  is  the  phase  velocity 
of  wave  with  angular  frequency  O'  ,  S  is  the  angle 
between  wind  and  waves  and  Uc  is  the  convection 
velocity  of  the  turbulent  eddies  with  wave  number  k. 
This  velocity  is  generally  considered  to  be  about  equal 
to  the  mean  wind  velocity  W  at  a  distance  2R/H  above 
the  sea  surface. 


The  Phillips  mechanism  calls  for  wave  energy  to 
grow  linearly  with  time  and  should  be  particularly 
important  for  waves  moving  with  approximately  the  same 
speed  and  direction  as  the  wind.  The  coefficient  of 
this  linear  growth  will  be  denoted  by  the  parameter  ci  » 
oL(f,  9  ,x,y,t}.  It  has  been  shown  by  Hasselmann  (1960) 
that 


P  (  k.  O) 


(2.3c  > 


26 


wher*  P(k,(r)  is  the  three  dimensional  spectrum  of  the 
random  atmospheric  pressure  fluctuation,  Fm  is  the 
density  of  sea  water,  g  is  gravitational  acceleration 
and  has  unit  of  (length/.  Barnett  (1968)  has 
proposed  a  suitable  function  for  Pdc,^*)  based  on  the 
work  of  Priestley  (1965).  The  procedures  are  as 
follows . 

Let  P(x,t)  be  the  fluctuating  static  pressure 
at  time  t  auid  point  (x,y).  Define  the  co-variance 
function 

R  Cr,  r)  =  <  ?(x+r,  ttt)>  IXJD 


where  the  brackets  denote  an  ensemble  average. 
The  quantity  desired  is 


=  <2K)'^]  I  f?(r.C)Cos(K.f+>L)<ifciL 
'c 


(  2.32) 


Priestley  defines  the  cross  spectral  density  as 


B  (r.a-) 


=  f 


I?  (  r  .C  >  exp  (-i(>i)  ol  L 


=  Co  (  r.o  )  +  t.  2  (  r.  b*)  (2.33) 
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whsra  c,  and  q  are  referred  to  as  the  Co-end 
Ouedreture-apectruRi,  respectively.  Priestley 

eventually  arrives  at  emprical  representations  for  c* 
and  q  as 

C*  ir.o.)  «  ^co.)  ejfp  I  -  I  r,  (  -  lr,|]  at, 

{2.3^) 

and 

9.  Cr.O")  **  C>)  I  -  V'l  ir,  t  -  I  Tal  j  3*n  aT, 

C2.i5) 

where 

A  ■  (1/tn) 

W  is  wind  speed  (m/sec) 
y'l  ■  0. 33  A*  **  (i/ra) 

-  0.52  (l/ra) 

<HP‘)  «  5.13  X  (0'‘*w‘/>* 

is  pressure  power  spectrum  (m  -sec) 

r,  and  r*  are  distance  associated  with  unit 

vectors 

% 

r,  and  ra  directed  downwind  and  cross wind, . 
res pectively . 

Substituting  Equation  (2.33),  (2.34)  and  (2.35)  into 
(2.32)  and  choosing  the  real  part  of  P,  then 

I  exp  I  -Pilr,l  j  coi  Ck.r-Ar,)cjr 

C2.3t) 


P(K.  O*)  * 


I 
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Dicact  integration  yields 


P  (k.o-) 


r  1 

(  1 

i.2.37) 


where  S  is  angle  between  wind  and  wave  vectors . 
Barnett  also  compared  the  oC  value  with  experimental 
data  which  fit  reasonably  wall. 


This  resonance  mechanism,  as  the  name  implies , 
should  be  of  particular  importance  for  waves  moving 
approximately  in  the  same  direction  with  the  same  speed 
as  the  wind  field.  We  have  found,  as  will  be  shown 
later  that  while  this  mechanism  is  important  in  the 
initial  stage  of  wave  generation  and  also  contributes  a 
large  portion  of  the  energy  in  the  eventual  spectral 
peak.  It  cannot  account  for  the  energy  growth  for  the 
rest  of  the  spectrum. 


2.4. 1.2  Miles'  Instability  Mechanism 

Miles  (1957,  1959a, 1959b, 1962 )  explains  the 
energy  transfer  from  the  wind  to  the  waves  by 
hypothesizing  a  coupling  between  the  existing  waves  and 
their  induced  surface  pressure  fluctuations.  The  waves 
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grow  exponentially  in  time  with  the  spectral  peak 
shifting  toward  lower  frequency.  Of  critical 
importance  in  this  particular  theory  is  the  curvature 
of  the  wind  profile  at  a  height  above  the  sea  surface 
where  W  >  C.  The  coefficient  of  this  exponential 
growth  is  denoted  by  the  parameter  P  (f,9,x,y,t). 
Unfortunately,  the  theoretical  prediction  of  Miles ' 
theory  did  not  agree  with  the  measurements.  Snyder  and 
Cox  (1966)  recommended  a  modified  expression  of  p. 


S  -  s  C  k  •  w  -  0-) 


where  s  is  the  ratio  of  the  densities  of  air  and  water, 
and  is  the  wind  velocity  measured  one  wave  length 
above  the  mean  sea  surface.  Barnett  (1963)  rec«nmended 
a  slight  change  in  the  ratio  to  yield 


0  .  5sf  (  W  C 


where  the  wind  W  is  considered  19.5  m  above  the  sea 

surface.  One  can  make  an  argument  to  the  effect  that 

most  instability  theories  require  waves  moving  faster 

u/ 

than  the  wind,  L.e.,  •^<1,  either  to  experience  no 
growth  or  to  be  attenuated.  The  dot  product 


3i 


t 
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representation  in  (2.33)  and  subsequent  form  of  (2.39) 

are  just  simple  math«natical  expressions  to  fulfill 

this  requirement.  However,  after  the  waves  have  moved 

away  from  the  generating  area,  the  actual  form  of  for 

w 

•^<1  then  becomes  important.  If  ^  is  relatively  large 
and  negative,  substantial  dissipation  would  occur 

continuously  outside  the  generating  region.  A  similar 
circumstance  would  occur  under  head  wind  conditions, 
(i.e.,  7C/2  <  S<  2K/2) ,  on  the  other  hand,  if  ^  is  small 
but  generally  positive,  as  has  been  suggestcKi  by 

Phillips  (1966),  then  slow  but  continual  growth  would 
be  expected  in  the  presence  of  light,  favorable  winds . 
To  make  a  choice  between  the  two  possibilities,  it 
would  be  necessary  to  observe  the  behavior  of  spectral 
components  over  great  distance  or  time.  The 
measur^nents  of  Snodgrass  et  al.  (1966)  indicated  that 

for  long  waves,  at  least,  the  deuaping  is  negligible. 

It  is  to  be  assumed  that  3*0  Wcos  ( -^)  <0. 90. 

Both  wave  growth  mechanisms  suggested  above  are 
for  deepwater  waves .  It  is  not  clear  whether  they  are 
applicable  for  shallow  water  case.  In  the  present 
model  the  conbined  Phillips'  and  Miles'  mechanism  is 
used  anyhow  for  lack  of  a  better  model,  that  is, 

e,  +  $,-</+  (3  A 
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To  illustrate  these  two  mechanisms,  a 
simple  case  is  tested  here  for  a  steady  state,  one 
dimensional,  sero  current,  and  frictionless  situation. 
Equation  (2.29)  then  becomes 

Cjc,3© -jL  A(f,e)  -  o<  +  3ACf,e) 

Equation  (2.40)  is  a  linear  boundary  value  problem. 
The  solution  can  be  easily  obtained  with  given  boundary 
condition.  Here,  the  boundary  condition  is  set  as 
A(f,9)*0  at  x«0.  The  solution  for  the  spectral 
component  along  the  x^axis ,  A(f,0)  is  simply 

A(f.o)  - (exp -i]  a4i) 

Three  different  conditions  are  illustrated  here. 

1)  Onshore  wind  blows  over  a  sloping  bottom  at  right 
angle  to  the  shore  line.  First  of  all,  the  growth  of  a 
specific  wave  component  (T^Qsec)  with  respect  to  fetch 
is  examined.  Here  as  shown  in  Figure  2.2(a),  the 
growth  of  each  component  will  undergo  three  stages.  In 
the  initial  stage,  wave  will  grow  linearly  with  respect 
to  distance  and  then  grow  exponentially.  Eventually, 
the  wave  will  reach  the  equilbrium  range.  Figure 
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2.2(b)  illustratas  the  fully  davelopad  spactrum  in 
shoaling  watar  at  location  “I**  basad  upon  tha  combined 
Phillips'  and  Milas'  mechanism.  The  major  peak  of  tha 
spectrum  shifts  tovrard  tha  wind  frequency 
(fM*g/2XW).  After  reaching  it  then  shifts  slowly 
toward  even  lower  frequency.  The  light  bump  near 
f«0.02  HZ  is  due  to  the  discontinuous  'cut  off  of  the 
Miles  instability  mechanism.  While  such  a  sharp 
cut-off  is  unrealistic,  it  will  not  significantly 
affect  subsequent  results.  The  fully  developed  P-M 
spectrum  for  the  same  wind  condition  is  also  shown  in 
Figure  2.2(b)  for  comparison  purposes.  The  Phillips' 
and  Miles '  mechanisms  allow  for  a  certain  amount  of 
growth  in  long  wave  components  which  do  not  appear  in 
the  fully  developed  P-M  spectrum. 

2)  offshore  wind  blows  over  a  sloping  bottom  as  shown 
in  Figure  2.3.  The  growth  of  a  specific  wave  component 
with  respect  to  fetch  still  maintains  the  three  growing 
stages  but  the  full  spectrum  at  the  location  "I"  is 
quite  different  from  the  above  case  for  lack  of 
shoaling  effect.  The  comparison  between  the  fully 
developed  P-M  spectrum  and  the  spectrum  generated  by 
Phillips'  and  Miles'  actually  have  very  similar  shape. 
The  peak-energy  frequencies  also  coincide  with  each 


oth«r.  The  total  energy,  thus  the  wave  height,  is 
apparently  lesser  for  the  shallow  water  case. 

3)  Wind  blows  over  a  deep  uniform  bottom.  The 
spectrum,  following  the  same  pattern  as  case  (2),  is 
shown  in  Figure  2.4.  Although,  at  present,  there  is  no 
field  verification  on  the  applicability  of  Phillips' 
and/or  Miles'  mechanism  in  shallow  water,  the  results 
given  above  do  appear  reasonable. 

2.4.2  Wave  Dissipation 

The  known  wave  dissipation  mechanisms  are 
bottom  friction,  wave  breaking,  percolation,  and  wave 
induced  bottom  motion.  The  former  two  are  usually 
considered  to  be  the  doninating  factors  and  both  are 
non-linear  mechanisms.  The  latter  two  proesses  are 
minor  and  linear  mechanisms.  In  this  reseeurch,  only 
the  two  major  mechanism:  bottom  friction  and  wave 
breaking,  and  depth  limited  wave  are  considered.  The 
two  major  mechanisms  are  the  term  in  the  energy 
transport  equation  whereas  the  last  term  limits  the 
growth  of  the  total  energy  under  the  spectral  curve  in 
shallow  water. 


2. 4. 2.1  Bottom  Friction 


37 


Energy  dissipation  due  to  bottom  friction  is 
primarily  the  result  of  two  physical  mechanisms:  one 
involves  the  dissipative  work  done  against  turbulent 
shear  stresses*  the  so-called  turbulent  bottom 
friction,  induced  by  water  particle  motions  near  the 
bottom  and  the  other  is  responsible  for  the  dissipative 
work  done  against  the  viscous  forces  induced  at  the 
permeable  bottom  where  percolation  occurs . 

In  general,  turbulent  bottom  friction  is 
dominant  over  other  dissipative  processes  when  the 
sediment  is  composed  of  sand  with  mean  diameter  in  the 
range  of  0. 1-0.4  tnn  for  which  low  permeability 
prohibits  percolation  and  viscous  dissipation  is 
negligible.  Putnam  and  Johnson  (1949)  investigated 
this  problem  first.  They  used  the  quadratic  friction 
law  to  derive  the  r^te  of  energy  dissipation  for 
sinusoidal  surface  wave.  Assuming  Gauss ian-distributed 
surface  wave  field,  Hasselmann  and  Collins  (1968) 
derived  the  rate  of  energy  dissipation  for  a  random 
sea.  Their  derivations  are  reviewed  briefly.  The 
bottom  friction  can  be  represented  by  the  quadratic 
law. 

*  -  Pw  Cf  luel  U.45) 
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where 

(i«l,2)  represent  the  two  horizontal 
components  of  water  particle  velocity 
at  the  bottom, 

Cfi  is  the  bottom  friction  coefficient. 

The  total  dissipation  of  the  wave  field  is  then  given 
by 

$  -  <7b;  U*.  >  —  -  P.*  Cf  <  llu|»>  (2.43) 


<  >  denotes  an  ensemble  average  over  all  possible  wave 
components  in  a  given  wave  spectrum.  The  repeated 
subscript  represents  summation  over  the  field  i«l,2. 
After  Hasselmann  and  Collins'  result  has  been  changed 
from  wave  number  space  to  directional  frequency  space, 
^  becomes 


$  (f.e) 


3  k* 
cos  h  *  k  h 


<u,>  A(f.e) 


€  A  (-f.e) 


(2.44) 


with 


<<.<►> 


K*  kh 


(2.4  5  ) 
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where  A(f)  is  energy  apectrum  given  by  Equation  (2.17). 
The  above  expressions  are  used  in  the  noncoupled  case 
because  the  basic  assumption  of  this  research  is  linear 
process  . 


As  the  bottom  friction  coefficient, 
Bretschneider  and  Reid  (1954)  evaluated  the  coefficient 
to  be  c^»0.01  for  the  sandy  bottoms  of  the  Gulf  of 
Mexico.  Hasselmann  and  Collins  (1963)  c^  to  be  0.015 
using  the  wave  spectra  measured  offshore  of  Panama 
City,  Florida.  Based  on  these  two  studies,  c^  values 
of  the  order  10  have  been  widely  used.  However,  the 
friction  coefficient  has  been  observed  to  vary 
significantly  above  and  below  this  range  depending  on 
the  sand  grain  dieuneter  and  existence  of  bottom 
ripples.  Laboratory  and  semi-theoretical  studies  by 
Jons  son  (1965)  indicate  that  Cf  is  a  function  of 
Reynold's  number  and,  also,  of  relative  roughness. 
Dingier  (1975)  and  Nielsen  (1977)  individually  found 
the  criteria  for  existence  of  bottom  ripples  and,  also, 
for  predicting  the  lengths  and  heights  of  such  ripples. 
Based  on  these  results,  a  method,  assuming  the 
roughness  height  being  equal  to  the  sediment  diameter, 
was  proposed  by  Hsiao  and  Shemdin  (1978)  to  estimate 
the  botttom  friction  coefficient  over  a  broad  range  of 
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conditons . 

2.4. 2. 2  Wave  Breaking 

It  is  a  well  observed  fact  that  under  sustained 
wind  action,  waves  will  grow  to  some  limiting 

amplitudes  then  breakings  or  whitecaps  will  limit  their 
growth.  This  phenomenon  can  be  handled  in  practice  by 
multiplying  the  wind  generation  term  ol  +|!J\  by  a 

weighting  factor  ju  where  M  is  unity  during  generation 
and  zero  at  breaking.  The  condition  Is  stated  as 

&,  +  &,  =  U  +  iSA) 

where  X(  is  1  for  A  unsaturated?  44  is  0  for  A 

saturated.  Phillips  (1958)  proposed  that  in  the 

saturated  range  of  wave  spectrum  in  deep  water  the 
expression  is 


A  a,  §),.w..i  =  B3V*@(,s)  u-*?; 

Barnett  (1968)  defined  (§)  as 

0  CS)  *  (2.48) 


i 


The  coefficient  ”  B  "  has  been  determined 


experimentally  by  Burling  (1959),  Hicks  (1960),  Kinsman 
(1960),  Pierson  et  al .  (1962),  Longuet-Higgins  et  al . 

(1963)  and  others,  as  a  value  varying  between  0.0080 
(Longuet-Higgins  et  al)  to  0.0148  (Burling)  with  an 
average  value  about  0.0123.  The  value  0.0123  is  used 
herein  (Phillips,  1966). 

In  shallow  water,  Kitaigorodskii  et  al .  (1975) 

have  extended  Phillips  arguments  to  case  of  finite 
depth  and  obtain  a  modified  equaion  of  the  following 
form: 

A  (-f.  5  )  =  B3V'‘©14)  r  (w„)  U.4<»> 

where  P  is  a  non-dimensional  function  of  the  quantity 
'^)^  .  The  function  varies  from  1  in  deep  water 

to  0  in  depth  h»0.  When  Wh  is  less  than  1, 

r tWh)  «  “  wh* 

Therefore 

A  5)  =  03hO.-^®  (S) 
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..Observational  evidence  for  this  form  has  been 
reported  by  Thornton  (1977),  Ou  (1980),  and  the  Vincent 
(1981)  in  the  laboratory  and  the  field. 


Finally,  when  waves  reach  shallow  water  they 
will  break  because  of  depth  limitation.  The  energy  is 
supposedly  saturated  but  Eq  (2.50)  can  not  account  for 
all  of  them  (the  low  frequency  components)  since  it 
opens  at  low  frequency.  In  the  present  model,  an 
empirical  criteria  based  upon  total  wave  energy  in 
terms  of  significant  wave  height  is  employed  such  that 
(Divoky  at  al . ,  1970) 


C 


b 


O  .  I  2  tan  K  C  k  fi )  b 


where  H3*4Je  with  E  the  total  wave  energy  under  the 


spectrum. 


CHAPTER  3 


NUMERICAL  SCHEME  OF 

STATIONARY  MAVE  SPECTRAL  TRANSFORMATION 


After  winds  steadily  generate  waves  for  quite  a 
long  period,  stationary  wave  spectral  transformation 
will  be  reached.  Without  knowing  the  processes,  the 
shallow  water  wave  spectrum,  at  any  specified  location, 
is  then  simply  the  contribution  resulting  from  energy 
associated  with  different  frequency  bands  of  the  deep 
water  wave  spectrum. 

Based  on  the  energy  flux  conservation,  Shiau 
and  Wang  (1977)  developed  numerical  model  of  stationary 
wave  spectrum  over  irregular  bottom.  The  present 
study,  based  on  the  same  numerical  schone  as  Shiau  and 
Wang,  add  wave  generation  and  wave  dissipation  terms  in 
the  stationary  wave  spectral  transformation  to  fit  the 
real  phenomenon  in  the  ocean  environments . 

3.1  Basic  Equations  for  Stationary  Wave  Spectral 
Trans  formation 
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Under  the  assumption  of  steady  state,  basic 
equations  for  wave  spectral  transformation,  i.e.. 
Equations  (2.21),  (2.25),  (2.28)  and  (2.29)  become 

i)  Wave  Conservation  Equation: 

^  =*  O 

or 

»  cortstort-i  or  C3. 1) 

ii)  Wave  Number  Conservation  Equation: 

t  9ktartli  kh  )^+  UkcosS  +  \/kSi*n  9  =»  -4^  (3.2) 

iii)  Wave  Number  Irrotational  Equation: 
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di.  is  Phillips  resonant  mechanism 
(3  is  Miles  instability  mechanism 
M  is  breaking  coefficient 
{  is  bottom  dissipation 


Eq  (2.30) 
Eq  (2.39) 

Eq  (2.44) 


3.2  Numerical  Procedure 

Based  on  Equations  (3.1),  (3.2),  (3.3)  and 
(3.4),  Noda,  et  al.(1974)  developed  numerical 
procedures  to  compute  the  wave  angle  and  wave  height  as 
functions  of  space  variables  for  beaches  with  periodic 
bottom  variations  in  the  longshore  direction.  The 
advantage  of  Noda's  method  is  that  it  predicts  the  wave 
angles  and  wave  heights  at  certain  points  rather  than 
along  a  wave  ray.  This  procedure  leads  itself  well  to 
use  in  the  finite  difference  model  because  calculations 
are  performed  at  points  which  lie  in  the  center  of 
rectangular  grid  elements  which  are  part  of  a  larger 
grid  mesh.  The  basic  idea  to  get  a  steady  solution 
converged  by  finite  difference  method  is  to  use  the 
relaxation  scheme  which  will  be  shown  in  detail.  It  is 
necessary  to  discuss  the  error  estimates  before  the 
numerical  scheme  to  be  used. 
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The  acccuracy  of  the  finite  difference  method 
is  mainly  limited  by  three  factors .  First  there  is  the 
truncation  error  which  is  the  difference  between  the 
true  solution  of  the  original  equations  and  the 
approximate  solution  of  the  finite  difference  equation 
which  replaces  the  differential  terms  by 
finite-difference  terms.  The  truncation  error  can  be 
eliminated  as  the  finite  difference  mesh  size 
approaches  zero.  Secondly,  there  is  the  round-off 
error  caused  by  the  limited  number  of  digits  stored  by 
the  computer  for  each  number.  As  the  computer  used 
includes  12  significant  digits,  no  appreciable  errors 
were  expected  from  this  source.  Thirdly,  there  is  the 
residual  error  caused  by  the  differen.ee  between  the 
exact  solution  to  the  finite  difference  equations  and 
the  actual  solution  reached  when  iteration  scheme  is 
stopped.  Obviously  if  a  large  number  of  iteration 
cycles  were  used  this  error  could  be  reduced  virtually 
to  zero  but  economy  insists  that  the  program  be  stopped 
as  soon  as  the  residual  error  is  reasonably  small. 

To  facilitate  numerical  computation,  the 
governing  equations  are,  out  of  necessity,  to  be 
transformed  into  finite  difference  form.  A  mixed 
forward,  backward  and  central  difference  scheme  is  used 


to  minlniz*  the  secondary  restrictions  at  the  expense 
of  nonuniformity  of  the  order  of  errors.  The  reference 
to  the  grid  schemes  is  shown  in  Figure  3.1. 

3.2.1  Wave  Kinematics  and  Boundary  Conditions 

The  wave  number  and  wave  angle  are  solved 
through  numerical  iterations  of  wave  number 
conservation  equation  and  wave  number  irrotational 
equation.  In  other  words,  for  given  spatial 
distribution  of  U,  V  and  h.  Equations  (3.2)  and  (3.3) 
enable  us  to  compute  k  and  Q  in  the  wave  region 
provided  the  boundary  condition  is  specified.  The  wave 
number  conservation  equation  becomes 

(SK;.a  iaeh  h.’.j)]  ^  +  ^-^5*4  •<'.4  ®  + 

V;.j  Kj.j  5.via%j  =*  -Ip  U-i") 

which  is  solved  through  a  Newtonian  iterative  technique 
defined  as 


k 


k  oJ  J 


■f  (kalj) 


13.4) 


where  f'(k)  is  the  first  derivative  of  f(k)  which  is 
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defined 


f  (k)  ■■  9ktarthkh  -  ^  -  (JkcojG  -  VkSinS  j  (3,7) 

The  iteration  is  considered  to  satisfy  the  relaxation 
criteria 


j  kftew  “■  kold  1  ^ 

O.oo  1 

1  k 

new  1 

(3.8) 

The  value  k  is  then 

fed 

into 

the 

wave 

nxnnber 

irrotational  equation 

to 

solve  the 

wave 

angle. 

Expanding  Equation  (3.3) 

yields 

^  -f  de  1  dk'l 

dx  ’  k  dyj 

+  0  ^ 

^  +  _L  .iil'l 
dy  K  dxj 

*0 

U.9) 

From  Equation  ( 3 . 2 } « 

dr 

dk 
'  dX 

can  be  derived  as 

dx 


“  {  k-^-  (UiinG  -  VeeSQ)  -  ^  (cose  "I—  +«J"9|7) 

cose  +  Vs.'n® 


-  9k^seeh^kh  dh  1  J.  [ 

2(3kto«hkW)^  dX  J  •  1  ^ 


3(ta»ihkh  4  Khsech^kh)  | 

2  ^Sktoohkh  )  ^  J 


U.lo) 


dk 

dy 


( 


k  Ci;5.v«0  -  vcose)  -  k  (cosG 


Jk 
2 


k*sech^kh  j.  I  ^  ^5.v,9 

^3kt«f.hkh)  ■'i  dy  J  *  I 

H>  kh  s»ch*  kh  )  1 

Wkh) ^  J 


3  (ioi»hkH 

2  ^  3k  i«« 


C3.  II) 
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However «  Equation  (3.9)  really  needs  -L 

K  w  if  K  iy 

which  can  be  represented  as  follows: 


I  r  Us:nd  -  Vc»s9  )  a  JL 

k  JX  *  JJ'  ^  AA  ^  ^x 

{3.12) 


\  ate  ae  i  Us.oO  >  v<^s_e  )  ^  j.  _i*L_ 

T  TT  *  ^  aa  '  K  dy 

{3.13) 


with 


AA  =  u«»9  +  Vj.os  +  i  (>■* 


^  -  U^®»6  “  Vs*"©  ) 


13.14) 


and 


i.  — 
k  dX  ” 


-  UK  case  -  VK  r.*^) 
S.’nh  t2Kh) _ 


AA 


Ci,«5) 


I  _£k 
K  dy 


((h-UkcgsG  -VK3.h€)  Jh 

~  77»»»>Q)  ~ _ s.nhcakh) 

AA 


C3./(>) 
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Substituting  (3.12)  and  (3.13)  into  (3.9),  the 
following  final  equation  is  arrived  at 


^0  ^  ^  S.'e©  (u  *.*•©- VC«5e) 'I  . 


dk 


s;n  6 


JL  dk 
k  *y 


Cos  G 


0 


O»9(0«»'«9  »Vc<t^)l 

AA 


j 


y.iij 


Finite  differencing  the  above  equation  by  using  a 
forward  difference  in  x  and  a  backward  difference  in  y, 
and  solving  for  6:,j 


Gu  = 


( 

f  JL 

9k 

8j.j 

9/ 

f 

S.*n©; 

«y  ' 

©;^.i 

^cos6 

••J 

Cos  ©.-.3  “  7  Tx  ■*' 


Cos  9;.  j 


A  A  ,'.j 


^  iUs.v.©;,]  -VC4>S©.-^]  - 


U.I8) 


where 


cos  Q.-,; 


A  A  ivj 


(Us*"d;.i  -  Vcos9:.^) 


Cos9;.i 


5«  ^  9 

AX 


52 


AA;,j  is  defined  in  equation  (3.14). 

—  ^  and  -J-  are  defined  in  Equations  (3.15)  and 

K  av  tc 

(3.16).  The  computation  is  carried  out  to  satisfy  the 
relaxation  criteria 


&  new/ 


ii.2o) 


The  boundary  condition  of  0  is  specified  at  the 
offshore  boundary  at  the  x-axis .  The  offshore  boundary 
value  is  the  input.  The  condition  of  the  x-axis  could 
be  specified  in  such  a  manner  that©;,,  *  0;,a  without 
introducing  significant  error  if  the  longshore  bottom 
variation  in  the  vicinity  of  the  x-axis  is  slow. 


will  not  converge  when  approaches 
zero.  Using  the  null-current  condition  as  a  guide,  one 
arrives  at  the  following  condition; 


■  ~  '  >  'teun  0m 


where  0m  is  the  expected  maximum  wave  angle.  For 

example,  if  the  maximum  wave  angle  is  expected  to  be 

- 

less  than  60*,  the  ratio  of  should  be  larger  than 

1.7. 


Th«  currant  and  depth  gradients  required  in  the 


computation  are  determined  by  central  difference,  i.e.. 


(3.?2) 


d  X, 


2  AX. 


where  ^  represents  U,  V  or  h  and  x„  is  either  x  or  y. 


3.2.2  Wave  Dynamics  and  Boundary  Condition 


The  energy  transport  Equation  (3.4)  can  be 
derived  as  the  following  form. 


or 


where 
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! 


^y  ~  ■*■  “5y“ 

Finite  differencing  the  above  equation  by  using  a 
forward  difference  in  x  and  a  backward  difference  in  y, 
it  becomes 


{u  * C3c.se)  t  (vtcji.ne) 


+  Aj.j  C  "  Q  ;.j  ~  ^  (3  ;,j  “  ^  i.i  ) ^  “  0 

I3.f4) 

Rearranging  Equation  (3.24) <  A;,j  can  be  represneted  as 


A.i 


(vtC3s;.e)  (ufCjc.e)  ■ 


ay 


C3.25) 

Here  all  the  values  of  U,  V,  and  9  are  evaluated  at 
(i. j) . 


Again,  Equation  (3.25)  is  solved  through  a  row  by  row 
relaxation  until 


I  Anew  -  Aold  I  ^  0*001  I  Anew  I 


(3.26  ) 


The  boundary  condition  is  also  defined  in  a 
similar  manner  in  that  the  shoreward  boundary  is 
specified  as  the  input  and  the  values  on  the  x-axis  is 
taken  the  same  as  the  next  column  (A-,,,  »  A,-,^  ). 
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On  the  shoreward  direction,  the  wave  will 
eventually  break.  Under  such  conditions,  the  wave 
breaking  dissipation  is  more  important  than  bottom 
friction  dissipation.  In  the  present  model.  Equation 
(2.51)  will  be  chosen  as  the  criterion. 

3 . 3  Examples 

A  plane  beach  of  parallel  bottom  contours,  with 
the  distance  of  25,000  m  between  the  offshore  boundary 
and  the  shoreline,  is  shown  in  Figure  3.2.  The 
coordinates  are  chosen  as  x-axis  toward  offshore, 
y-axis  alongshore,  and  z-axis  upward.  To  apply  the 
numerical  scheme,  along  the  positive  x-axis,  the  domain 
is  divided  into  five  sections,  with  Ax*5000  m,  denoted 
as  location  1  to  6,  and  along  the  positive  y-axis,  it 
is  divided  into  two  sections,  with  Ay»5000  m,  denoted 
as  location  I  to  III.  Three  pressure  gauges  are 
installed  at  location  6II  (boundary),  location  5II,  and 
location  2II  of  which  the  water  depth  are  25  m,  20  m, 
and  2  m,  respectively. 

There  are  three  examples  to  be  discussed:  i) 
energy  flux  conservation  case,  ii)  bottom  dissipation 
case,  iii)  local  wind  generation  case. 
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3.3.1  Energy  Flux  Conservation  Case 

The  energy  flux  conservation  case  considers  no 

energy  generation  or  dissipation,  i.e.,  2G:*0,  in  the 

Ul 

process. of  the  wave  spectral  transformation.  The  wave 
spectra  at  different  locations  are  shown  in  Figures  3.3 
and  3.4.  Since  the  shoaling  coefficient  K^h  • 
(unaffected  by  refraction), 

wave  height  in  shallow  water 
= - - - 

wave  height  in  deep  water 

I, 

is  less  than  1  when  0.06  *  0*5  and  larger  than  1 

when  <  0.06, (Lo  is  the  wave  length  in  deep  water), 

the  wave  height  will  increase  or  decrease  in  shoaling 
water  depending  on  the  range  of  .  In  the  example 

given  in  Figure  3.3,  the  energy  spectrum  at  location 
5II  is  lower  than  the  input  spectrum  (location  6II), 
because  the  value  of  of  all  the  frequency 

components  fall  between  0.06  and  0.5.  On  the  other 
hand,  the  components  of  energy  spectrum  at  location  2II 
are  lower  than  the  inputs  at  low  frequency  but  higher 
than  the  inputs  at  high  frequency.  Figure  3.4  shows 
the  refraction  effects  in  shoaling  water.  The 
directional  sensitivity  is  further  examined  in  Figure 


Figure  3.4  Variation  of  Ehergy  Spectrum  in  Energy  Flux 

Conservation  Case  at  Boundary,  Location  5II  and 
Location  211 
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3.5.  It  plots  ths  psrcentaga  error  of  spectral  density 
of  each  frequency  component  for  various  approaching 
angles  that  deviate  from  a  designated  mean  wave  angle. 
It  can  be  seen  that  error  increases  with  increasing 
wave  period.  More  importantly,  the  approaching  angle 
plays  a  significant  role  in  directional  sensitivity. 
For  waves  with  small  approaching  angle  (perpendicular 
to  the  shoreline) ,  the  wave  transformation  is  far  less 
sensitive  to  the  wave  angle  variations  than  for  waves 
with  large  approaching  angle  (parallel  to  the 
shoreline).  Therefore,  the  energy  spectrum  is 
dominated  by  frequency  variations  for  small  approaching 
angle  but  by  wave  angle  variations  for  large 
approaching  angle. 

3.3.2  Bottom  Dissipation  Case 

The  bottom  dissipation  case  concerns  the  energy 
loss  due  to  the  bottom  friction  in  the  process  of  wave 
spectral  transformation.  As  shown  in  Figure  3.6,  the 
influence  of  bottom  friction  is  quite  limited  at 
location  5II  because  few  wave  components  are  affected 
by  the  bottom.  The  difference  between  the  wave 
spectrum  with  bottom  friction  and  the  wave  spectrum 
without  bottom  friction  is  only  a  few  percent  for  the 


Figure  3.5 


Figure  3.6  Comparison  between  Energy  Spectrum  with  and 
without  Bottom  Friction  at  Location  51 I 


low  frequencies  under  consideration  and  these  two 
spectra  practically  overlap.  However,  in  Figure  3.7, 
the  influence  of  the  bottom  friction  is  quite 
significant  at  location  211  because  most  of  the  wave 
components  have  already  felt  the  bottom.  It  is 
reasonable  that  the  bottom  dissipation  effect  for  long 
wave  is  more  significant  than  that  for  short  wave  in 
shallow  water  because  long  wave  feels  the  bottom 
earlier  than  short  wave. 

3.3.3  Local  Wind  Generation  Case 

The  local-  wind  generation  case  deals  with  the 
energy  generation  in  the  process  of  wave  spectral 
transformation.  Considering  no  energy  input  at  the 
boundary,  the  steady  local  wind  generates  waves  in  a 
closed  region.  As  shown  in  Figure  3.8,  the  local  wind 
wave  spectrum  at  location  5II  still  stays  at  high 
frequency,  whereas,  with  increasing  fetch,  not  only  the 
local  wind  wave  spectrvim  at  location  2II  increases 
singificantly  but  its  spectral  peak  shifts  toward  the 
lower  frequency.  Hasselmann  et  al.  (1973)  in  the 
JONSWAP  measurements  which  traced  the  development  of 
the  wave  spectrum  with  increasing  fetch  under 
conditions  of  steady  offshore  winds  illustrate  this 


Figure  3.7  CompdcLson  between  Energy  Spectrum  with  and 
without  Bottom  Friction  at  Location  2II 
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phenomenon.  Since  the  Kitaigorodskii ' s  shallow  water 
equilibrium  energy  criteria  .open  at  one  end,  can  only 
handle  the  wave  breaking  for  a  specific  range,  the 
criteria  for  the  breaking  of  the  I'^ng  wave  components 
needs  to.  be  restricted  by  the  depth-limited  wave 
condition.  The  comparison  of  the  local  wind  wave 
spectrum  at  location  5II  and  location  2II  with  P-M 
spectrum  is  shown  in  Figure  3.9.  The  result  at 
location  5II  is  fair  because  the  energy  spectrum  of 
that  particular  location  hasn't  been  fully  developed 
yet.  The  result  at  location  2II  is  also  fair  because 
the  P-M  spectrum  is  only  valid  in  deep  water. 


Figure  3.9  Comparison  of  Energy  Spectrum  of  Local 
Wind  at  Location  511  and  21 I  with 
Pierson-Moskowitz  Spectrum 


CHAPTER  4 


NUMERICAL  SCHEME  OF 

NCN-STATIONARY  WVE  SPECTRAL  TRANSFORMATION 


In  the  previous  chapter,  the  numerical  scheme 
for  the  stationary  wave  spectral  transformation  was 
developed.  For  the  case  of  non-stationary  wave 
spectral  transformation,  the  problem  is  more 
complicated  because  the  kinematic  and  dynamic 
characteristics  change  with  time.  This  chapter  extends 
the  model  to  unsteady  state  case  when  the  offshore 
boundary  condition,  and  the  local  wind  effect  in  the 
domain  of  consideration  are  no  longer  stationary.  This 
is  a  common  situation  when  a  hurricane  or  storm  sweeps 
through  offshore  of  a  coastal  region. 

4.1  Basic  Equations  for  Non-Stationary  Wave  Spectral 
Transformation 

Basic  equations  for  non-stationary  wave 
spectral  transformation  have  been  developed  in  chapter 
2.  They  are  repeated  here; 
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$  13  bottom  dissipation  Eq.  ’(2.44) 

4.2  Method  of  Solution 

The  numerical  scheme  described  in  the  previous 
chapter  can  not  be  easily  extended  to  the  non-s tatioary 
case  because  the  nature  of  the  differential  equations 
changed  s ignif iciantly  with  the  addition  of  the 
unsteady  term  .  An  entirely  new  numerical  scheme 
has  to  be  developed. 

4.2.1  Two-Step  Lax-Wendroff  Scheme 

The  finite  difference  scheme  for 

initial-boundary  value  problems  can  be  classified  into 
two  categories :  one  is  implicit  method  and  the  other 
is  explicit  method.  Although  the  implicit  method 
ensures  unconditional  stability  criteria,  it  usually 
involves  lengthy  computations  due  to  matrix  inverse. 
Therefore,  an  explicit  method  is  explored  here.  The 
two-step  Lax-Wendroff  scheme,  i.e.,  the  first  half  step 
is  Lax-Fr iedr ichs  scheme  and  the  second  half  step  is 
leap-frog  scheme,  was  used.  The  reasons  for  selecting 
this  scheme  are  that  (1)  it  is  second  order  accurate, 
and  relatively  low  in  dissipation:  (2)  it  is 
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single-lavelad  and  eaay  to  apply;  and  (3)  it  I3  widely 
used  by  fluid  dynamicists  and  hence  it  would  be 
desirable  to  understand  fully. 

Before  applying  the  actual  problem  posed  above, 
the  numerical  scheme  is  explained  first  with  the 
application  to  an  initial-boundary  value  problem  which 
has  the  similar  characteristics  as  the  actual  problems 
but  with  exact  solution.  In  this  way,  the  validity  of 
the  numerical  scheme  can  be  established.  The  stability 
criteria  was  then  examined. 


The  example,  with  a  non-constant  coefficient 
and  non-homogeneous  term,  is  shown  as  follows . 

Governing  Equation  (  G.E.): 


Initial  Condition  (  I.C.  ); 

U  (x,0)  »  X  /  S9L62 
Boundary  Conditions  (  B.C  ): 

U(0,t)  =  0 


U(35,t)  -35  /  (S9L62  +  t) 
Then  the  exact  solution  is 


U(x,t)  -  X  /  (39162  +  f) 
Let's  define  the  coefficients  as 


C  ) 


Q  ,  (y,  t)  =  • 

Qi  IV,  t )  *  • 

Equation  (4.5)  becomes 


2y 

+  t 
y  -» 

4?‘42  +  t 


O,  {y,t) 


-  Qj  t  y,t)  U  i4.7) 


According  to  the  two-step  Lax-Wandroff  scheme,  the 
difference  equation  of  Equation  (4.7)  becones 


i)  Lax-Pr iedrichs '  Scheme: 

difference . 


modified  forward 


rt-i 


At 


n-l  n>l 

=  -  Qj  .•  u ; 


ii)  Leap-Frog  Scheme:  three  time  level. 


{4.8) 


H  +  l  rt-l  ,  \ 

U.'  -C;.,  ;  _  ^  .•» 

-  +  (3,  ;  - ;; -  *  -  Qi  (x.t-),- 

2^t 


{4.<R) 


where 


the  subscript  i,  j  represent  the  discr itization  of 
x-axls  and  y-axis, 

the  superscript  n  represents  the  discr itization  in 
time  axis . 

Using  a  time  step4t/2  instead  of  At  in  the  first  two 
half  steps,  the  accuracy  of  the  first  calculation  has 
been  improved.  The  computation  procedures  are 
illustrated  schematically  in  Figure  4.1  and  4.2. 

4.2.2  Stability  Analysis 

Richtmyer  and  Morton  (1967)  showed  that  the 
well-known  Lax  Equivalence  Theorem  is 

"For  a  well-posed  initial-value  problem  for 
linear  differential  equations,  and  for  a  consistent 
difference  approximation,  stability  is  necessary  and 
sufficient  for  convergence." 

In  other  words,  to  show  stability  analysis  is 
to  3  how  convergence . 


4. 2. 2.1  Homogeneous  Equation 


Figure  4.1 


Initial  Lax-Friedrichs  and  the  First  Leap-Frog 
Scheme  Solution  Steps 
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n+1 

n+2 

2at 

Sat: 

Lsap— Prog 


Leap-Frog  Solution  Scheme 
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The  homogeneous  equation  for  Equation  (4.7)  is 


+  Q.cx.t)  -4^  -  0 


l4.lo) 


The  stability  criteria  can  be  found  as 


|o.u.t)^l  4  I 


or 


Ai  4  I 


(X,  Cx,t) 


C4.i»  ) 


If  a,(x,t)  is  constant,  the  stability  criteria  will 
vary  with  Ax  only  and  follow  the  Equation  (4.11). 
Whereas,  a,  (x,t)  is  not  a  constant,  the  time  step  At 
actually  used  is  much  less  than  that  given  by  the  above 
criteria. 


4 . 2 . 2 . 2  Non-Homogeneous  Equation 

The  stability  criteria  for  the  non->homogeneous 
equation  has  to  meet  the  following  condition,  one 
condition  is  the  same  as  that  of  homogeneous  equation, 


i  .e . , 


Ai  1  i 


AX 

a.cr.t) 


77 


the  other  condition  is 

(Oa  CK.  t)  1  <  I 


(4./2) 


The  domain  of  the  problem,  with  distance  of  35  along 
the  x-axis,  is  divided  into  five  sections,  denoted  as 
location  1  to  6.  By  using  ^t  equal  to  84  and  Ax  equal 
to  7,  the  comparison  between  the  exact  and  the 
approximate  solution  from  location  2  to  5  are  shown  in 
Figure  4.3.  It  seems  that  the  approximate  solution 
becomes  unstable  after  a  few  time  steps.  The  reason  is 
due  to  the  time  level  of  the  non-homogeneous  term  in 
Leap-Frog  Scheme.  By  using  (n-l)th  time  level  instead 
cf  nth  time  level  for  the  non-homogeneous  term,  the 
problem  will  be  handled.  Equation  (4.9)  is  rewritten 
as  follows: 


ur-ur 


2At 


+  Q,U.i); 


*»  ( “  U;-t) 


2ax 


rt-l  A-l 

QaCX.t);  \J, 


(4.13) 


The  new  results  are  shown  in  Figure  4.4  and  agree 
reasonably  well.  Thus  the  effectiveness  v-)f  the 
numerical  scheme  for  solving  a  transient  problem  is 


demonstrated . 


4.2.3  Numerical  Dissipation  and  Dispersion  Analysis 


It  is  well  known  that  numerical  schemes 
introduce  numerical  dissipation  and  dispersion,  in 
roughly  the  same  way  as  physical  dissipation  and 
dispersion  occur  in  phenomena  of  fluid  flow. 
Dissipation  means  the  damping  of  Fourier  components  or 
plane  waves.  Dispersion  means  the  propagation  of  plane 
waves  of  different  components  at  different  speeds  c. 

For  the  case  of  one  dimensional  shallow  water 
wave  in  uniform  depth  and  zero  current,  the  wave 
conservation  equation  and  the  energy  transport  eqution 
can  be  written  as 


aA 


+  C 

+  C 


dK 

dX 


o 


o 


(4. 14  ) 


Both  equations  are  wave  equations  representing  a  wave 
of  form  k  or  A  propagating  without  change  of  shape  at 
speed  c  in  the  positive  x  direction.  To  illustrate  the 
effect  of  numerical  dissipation  and  dispersion,  an 
example,  with  a  function  f  representing  k  and  A  in 
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Equations  (4.14)  and  (4.15),  is  shown  as  follows:. 
Governing  Equations  (G.E.): 

_if.  +  c  -±L  -  o  u.ft) 

Initial  Condition  (I.C.): 

■f  Cv.c)  =  o 


Boundary  Condition  (B.C.): 


-f  ^o.t) 


o 
2t 
4- t 
O 


The  exact  solution  is 


t.<o 

O  <  t<^ 

2  i  t  <4 

t.  24 


-p  (.  -  C  -fe) 


which  is  recognized  as  having  the  wave  character  noted. 
According  to  the  two-step  Lax-Wendroff  scheme,  the 
differtsnce  equation  of  Equation  (4.15)  becomes 
i)  Lax-Friedrichs '  Scheme: 


At 


f... 


) 


+  c 


f. 


o  -4 
4» 


o 


(xf-,  IS  } 
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ii)  Leap-Frog  Schema: 


+  C 


2Ay 


O 


U./?) 


The  comparison  of  exact  and  numerical  solution  for 
r(»C"*r)»l/2  is  shown  in  Figure  4.5.  It  can  be  seen 
that  both  the  dissipation  and  dispersion  effects  are 
more  significant  for  At»l/2  than  for  At*l/4.  A  large 
number  of  spurious  oscillations  near  sharp  transition 
is  due  to  the  conservative  scheme  which  will  be 
discussed  later. 

To  illustrate  the  dissipation  and  dispersion 
affects  of  the  Lax-Fr iedr ichs  scheme  and  the  Leap-Frog 
scheme,  the  amplification  factors  for  the  two  schemes 
need  to  be  found.  By  using  the  Von  Neumann  method,  a 
funtion  can  be  represented  by  a  finite  Fouries  series. 


JL  ^ 

Z 

JZz.O, 


U 


i-4,  2^  ) 


where 

1  represents  the  wave  component, 
L  is  wave  length. 


-^.ja 


f  la.CI 


•<ja 
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Sub«tuting  Equation  (4.20)  into  Equation  (4.13), 
raapastivaly,  ona  can  obtain 


where  amplification  factors  for 

Lax-Priedricha '  scheme  and  leap-frog  scheme  become 
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Th«n  th«  amplltud«s  iRAl  for  th«  two  «sh«m««  ar« 
raspoctlvol/  by 

[f?Al  «  C 

and 

iRAl  -  I 

The  phase  angles  <p^  for  the  two  schemes  can  be 
respectively  as : 


“*  (  r  tous  ny) 


and 


C/)^  s  tool 


.1 


r$'*n 


2'KJ^ 


2K^ 


.X  ) 


Then  the  dispersion  effect  can  be  shown  as 

C  Vuw  __  ^A/ 

C  fLta  ^ 


given 


u.2t; 

shown 


^4,2^) 


wav«  length. 


represents  the  number  of.  samples  per 


The  dissipation  and  dispersion  for  the 
Lax-Friedrichs  scheme  and  leap-frog  scheme  are  shown  in 
Figures  4.6  and  4.7,  respectively.  It  can  be  seen  that 
both  dissipation  and  dispersion  effects  are  more 
significant  for  short  wave  than  long  wave.  The 
Laz-Friedrichs  scheme  is  dissipative  and  dispersive, 
%diile  the  leap-frog  scheme  is  non-dissipative,  or 
conservative  and  dispersive. 


4.3  Numerical  Analysis  for  Wave  Conservation  Equation 


Wave  conservation  equation  is  a  homogeneous 
equation  which  can  use  the  two-step  Lax-Wendroff 
scheme.  However,  it  can’t  be  solved  singly  because  it 
contains  angular  frequency  O',  wave  nvonber  1c ,  and  wave 
direction  9 .  The  wave  kinematic  equation  including 
wave  conservation  equation,  wave  number  conservation 
equation,  and  wave  number  irrotational  equation  have  to 
be  solved  simultaneously. 


4.3.1  Numerical  Procedures 


Figure  4.7  Dissipation  and  Dispersion  for  the 
Leap-Frog  Scheme 
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To  BOlv«  it  algobraically,  tha  vector  form  of 
wave  conservation  equation  has  to  be  decomposed  into 
two  scalar  forma,  one  in  x-direction  and  the  other  in 
y>dir action.  If  the  wave  direction  &  measured  from  the 
positive  x-axis  in  a  counterclockwise  direction. 
Equation  (4.1)  becomes 


x-direction 

d(Kc*se)  .  ^ 

U.3#) 

H  *  »x 

y-direction 

(4.30 

dt 

other  two  wave  kinematic  equations  remain 

to  be: 

(  9  kianh  kh  ) 

^  +  UKceS  e  Vk5.‘n0  - 

-p  (4,2^ 

and 


^(Kc0S9) 

^(KS>'n9)  _  ^ 

(43) 

For  given  spatial  distribution  of  U,  V,  and  h;  ^  ,  k, 
and  9  can  be  obtained  by  solving  Equations  (4.30), 
(4.31),  (4.2)  and  (4.3)  simultaneously  with  specified 
initial  and  boundary  conditions.  Using  the  two-step 
Lax-Wendroff  scheme. 


B 


1)  Lftx-fri«drlshs '  Sch«m« 
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Because  the  coefficient  is  a  constant,  the  time  step  At 
obeys  the  stability  criterion  as  given  by  Eq.  (4.20). 

4.3.3  Initial  Condition  and  Boundary  Conditions 

To  solve  the  wave  conservation  equation,  it  is 
necessary  to  choose  suitable  initial  condition  in  the 
domain  and  boundary  conditions  at  boundaries.  The 
initial  condition  can  be  decided  by  either  angular 
frequency  or  wave  number  and  wave  angle.  Based  on  the 
steady  state  assumption,  the  angular  frequency  does  not 
change  with  space.  It  is  constant  even  as  the  water 
depth  changes.  Hence,  the  initial  condition  is 
determined  more  conveniently  by  constant  angular 
frequency,  which  is  invariant  with  space,  than  by  the 
wave  nianber  and  wave  angle,  which  varies  with  space. 

As  for  the  boundary  conditions,  the  variation 
of  the  wave  energy  has  an  effect  on  the  offshore 
boundary  condition.  In  other  words,  the  kinematic 
condition  will  vary  if  the  dynamic  condition  changes. 
The  rest  of  the  boundary  conditions,  shore  boundary  and 
side  boundaries,  need  to  be  specified  to  carry  out  the 
numerical  scheme.  At  shore,  no  flow  of  water  occurs 
across  the  shoreline,  thus  the  shore  boundary  doesn't 
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vary  with  tha  changing  dynamic  condition.  The  side 
boundaries  require  that  the  longshore  bottom  variation 
In  the  vlnclnty  of  both  sides  Is  slow. 


4.4  Numerical  Analysis  for  Energy  Transport  Equation 

The  energy  transport  equation  is  a 
non-homogeneous  equation  which  is  the  same  type  as 
Equation  (4.7)  except  it  is  a  two  dimemslonal  equation. 


4.4.1  Numerical  Procedures 


After  expanding  and  rearranging.  Equation  (4.4) 
arrives  at 

4  Cu  +  Cj  ce$9)  4  tV+CjS.n©)  -|5  + 

A  (  v;  +  Cj  case)  4  w+qj*-"©)  4  I7  + 

4  (Xy,  ^  Oyy  ^  -M(3  -  €  ]  -  Xioi  * «  {4.37) 

or 


where 


^^4  C  V;4Cjc*j9)  15  4  ^V4CJ$.«9)  ^  =  C®  +MAt«‘)A  + 

U  i$) 

GL  *  tCjS.nO  -  cose 


-  Cj<ei0  -jp  -5.'e9  -3^  -j 


1 
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Applying  the  two-3t*p  Lax-Wendroff  9ch«m«,  bh«  finite 
difference  form  of  Equation  (4.33)  is 
i)  Lax-Friedrichs '  Scheme 


n  I  •*-»  •»-!  rt-t  \ 


4t. 


U  +  Cj  C0S6); 


1%  «l 


A  •! 


CA.-w. :  -  A  ) 

2^y 


UwH  rt«i  . 

•  **♦<  "*  A  I .  i-)  ^ 


2ay 


.  _  «»«l  r.-. 

vQ.-.;  +  ^p:.%  )  Ai*»i  +  A*»c^;.j 


**-•  .  A  •! 


•  •I 


ii)  Leap>-Frog  Scheme 
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To  avoid  instabilities,,  the  non-homogeneous  terra  has 
to  be  put  in  (n^-l)  time  level  as  shovm  in  Equation 
(4.40) . 


4.4.2  Stability  Analysis 

According  to  the  stability  analysis  for 
non- homogeneous  equation,  the  stability  criteria  for 
the  present  difference  equation  has  to  meet  the 
following  two  conditions: 

i)  . 
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In  general,  the  first  condition  can  be  stated  in  the 
following  manner:  the  speed  of  propagation  of  any 
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disturbance  in  the  model  must  be  less  than  or  equal  to 
the  speed  it  takes  the  disturbance  to  cross  one 
computational  grid  block  in  one  computational  time 
step.  The  disturbance  speed  is  in  general  the  speed  of 
gravity  wave  plus  the  time- independent  mean  current. 
Therefore,  in  shallow  water.  Equation  (4.41)  can  be 
expressed  approximately  as 


I  Ul  +  c 


In  general,  the  maximum  magnitude  of  the  wave 
speed  far  exeeds  the  current  speed,  the  stability 
criteria  can  be  expressed  as 


mejf 


which  is  a  type  of  two  dimensional  Courant  number. 

As  stated  before,  the  time  step  actually  used 
in  the  numerical  model  is  much  less  than  that  given  by 
the  criteria  above  due  to  the  fact  that  the  coefficient 
is  not  a  constant  in  Equation  (4.41). 
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4.4.3  Ini^tial  Condition  and  Boundary  Conditions 

It  is  necessary  to  specify  the  initial  and 
boundary  conditions  to  solve  energy  transport  equation. 
Commonly,  the  initial  condition  falls  in  one  of  the  two 
types :  one  is  for  calm  sea  where  the  wave  energy  is 
equal  to  zero  and  the  other  is  the  fully  developed  sea 
where  the  wave  energy  is  saturated. 

There  are  also  two  types  of  offshore  boundary 
conditions.  The  first  type  is  an  open  boundary.  It  is 
a  time  varying  boundary  condition  which  is  affected  by 
the  offshore  wave  condition.  In  general,  it  can  be 
prescribed  by  a  continuous  variation  of  wave  spectrum 
at  the  offshore  boundary.  The  second  type  is  a  closed 
boundary  which  is  not  affected  by  the  offshore  wave 
condition  but  self-adjusting.  In  order  to  apply  the 
numerical  scheme,  the  shore  boundary  condition  and  the 
side  boundary  conditions  also  need  to  be  specified.  At 
the  shoreline,  the  wave  numbers  are  to  be  perpendicular 
to  the  shoreline  and  the  bottom  slope  to  be  finite. 
The  side  boundaries  are  such  that  the  bottom  contours 
are  parallel  outside  the  region  of  computation. 


4.5  Application  in  Monochromatic  Wave 
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To  examine  the  mlcro-mechanlam  of  wave  spectral 
transformation,  the  cases  of  a  one-dimensional  and 
two-dimensional  monochromatic  vrave  are  illustrated 
here.  The  random  wave  model  and  the  field  data 
comparison  are  dealt  with  in  the  next  chapter. 

In  computing  wave  kinmnatics ,  the  initial 
condition  is  set  at  constant  angular  frequency  over  the 
entire  region  of  computation.  Therefore,  the  boundary 
conditions  should  also  be  specified  in  terms  of  angular 
frequency.  The  offshore  boundary  condition  can  be 
determined  from  a  continous  variation  of  wave  spectrum 
at  the  offshore  boundary.  The  shore  boundary  doesn't 
vary  with  changing  the  dynaunic  condition.  In  order  to 
proceed  with  the  numerical  scheme,  the  boundary 
conditions  at  both  sides  will  be  specified  as  follows: 
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In  the  energy  transport  equation,  the  steady 
state  wave  spectral  tran:> formation  is  chosen  as  the 
initial  condition.  Since  only  swell  condition  is 
considered  in  this  section,  an  open  offshore  boundary 
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condition  is  used.  The  boundary  condition  at  the 
shoreline  is  assumed  to  be  both  no  energy  flux  into  the 
shoreline  and  no  reflection  from  the  shoreline.  The 
boundary  conditions  at  both  sides  need  to  be  specified 
because  of  the  numerical  procedure.  They  can  be 
expressed  as  follows: 
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A  schematic  view  of  a  two-dimensional  case, 
with  horizontal  dimensions  60  m  by  150  m  and  water 
depth  1  m,  Ls  shown  in  Figure  4.8.  The  coordinates  are 
chosen  as  x-axis  toward  the  offshore,  y-axis 
alongshore,  and  z-axis  upward.  To  apply  the  numerical 
schone,  the  region  is  divided  into  five  sections  along 
the  positive  x-axls ,  with  Axs20  m,  denoted  as  location 
1  to  6,  and  it  is  divided  into  four  sections  along  the 
positive  y-axis,  with  Ay* I 5  ra,  denoted  as  location  I  to 
V.  The  time  step  At*l  sec  is  used  in  all  cases . 


4.5.1  One-Dimensional  Casa 


Schematic  View  of 
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In  this  ssction,  thr««  diff«r«nt  conditions 
will  bs  consldorsd.  Ths  first  cass  axaminss  th«  wav* 
transformation  process  whan  the  wave  at  the  offshore 
boundary  becomes  increasingly  shorter  (a  case  of 
diminishing  wind).  \  monochromatic  wave,  with  T*10  sec 
and  H«0.Lm,  propagates  along  the  negative  x-axis .  On 
reaching  the  steady  state,  the  frequency  of  the  wave 
component  stays  at  0.1  HZ  no  matter  where  the  location 
is.  As  shown  in  Figure  4.9,  the  three  plots  represent 
the  wave  spectrum  at  location  6III  (boundary),  location 
Sill  and  location  2III,  respectively^  the  %#ave 
spectrum  of  each  plot  denoted  by  t»0  sec  is  at  steady 
state.  The  wave  at  offshore  boundary  is  then  gradually 
adjusted  to  generate  a  higher  frequency  wave  component, 
which  increases  frequency  from  f«0.l0  HZ  to  f«0.l5  HZ, 
within  60  seconds.  When  t«20  sec,  the  frequency  of  the 
wave  component  at  location  61 II  is  higher  than  that  of 
the  wave  component  at  location  Sill.  However,  the  wave 
component  at  location  2III  is  still  at  steady  state. 
It  is  due  to  the  fact  that  the  component  at  the 
boundary  has  not  arrived  at  location  Sill  and  location 
2III  yet.  As  time  elapses,  this  phenomenon  occurs  also 
at  both  location  31 II  and  location  2III.  If  the  wave 
component  at  the  boundary  becomes  steady  again,  the 
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wav«  component  at  other  locations  continues  to  change 
until  it  reaches  the  new  steady  state. 

The  second  case,  shown  in  Figure  4.10, 
Investigates  the  vrave  transformation  for  a  shorter  wave 
with  T«2.S  sec  and  H«0.08  m.  Under  the  same  condition 
as  in  the  first  case,  all  the  short  wave  components 
will  follow  the  same  pattern  as  that  of  long  vrave 
component  except  that  the  short  wave  component  travels 
3lo%rar  than  the  long  wave  component  within  the  same 
time  period.  It  is  owing  to  the  fact  that  the  group 
velocity  of  a  short  wave  is  smaller  than  that  of  a  long 
wave.  The  short  wave  component  lags  behind  as  the  long 
wave  component  moves  ahead. 

In  contrast  to  the  first  two  cases,  the  third 
case  investigates  the  wave  transformation  when  the  wave 
at  offshore  boundary  becomes  longer  with  the  short  wave 
component  as  the  initial  steady  state  condition.  As 
soon  as  the  monochromatic  wave,  T«10  sec  and  Ha*0.1  m, 
propagating  along  the  negative  x-axis ,  reaches  steady 
state,  the  offshore  wave  is  adjusted  to  lower 
frequency,  from  f«0.10  HZ  to  f«0.05  HZ  within  60 
seconds.  This  phenomenon  is  shovm  in  Figure  4.11. 
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4.5.2  Two-Dimensional  Casa 

Three  different  cases  are  discussed  in  this 
section.  The  first  case  Illustrates  when  a 
monochromatic  wave,  with  T«10  sec  and  HaO.l  m, 
propagates  in  the  direction  of  9*20* with  respect  to 
negative  x-axls .  After  the  wave  transformation  reaches 
Its  steady  state,  the  wave  at  offshore  boundary  Is  then 
adjusted  to  generate  higher  frequency  component  from 
f*0.10  HZ  to  f«0.15  HZ  and  smaller  wave  angle  within  60 
seconds.  The  variation  of  energy  spectrum  In 
two-dlmenslons ,  as  3ho%im  In  Figure  4.12,  Is  similar  to 
that  In  one-dlmens Ion .  As  time  elapses,  as  shown  in 
Figure  4.13,  the  wave  component  approaching  the  center 
line,  i.e.,  6*0*,  from  right  to  left,  occurs  first  at 
boundary,  then  at  location  5III  and  finally  at  location 
2III.  It  is  because  the  wave  frequency,  wave  number 
and  wave  angle  are  correlated  during  the  non-s tatlonary 
process.  The  wave  angle  will  change  as  the  wave 
frequency  varies  with  time. 

Figure  4.14  provides  the  second  case  with  the 
same  condition  as  the  first  case  except  propagating  in 
the  direction  of  9»-20*  with  respect  to  the  negative 
X-axis.  As  time  elapses,  the  wave  component 
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approaching  the  center  line,  from  left  to  right,  is 
opposite  to  the  first  case. 

The  third  case  with  a  situation  just  opposite 
to  the  first  case  in  that  the  wave  at  offshore  boundary 
is  adjusted  to  generate  lower  frequency  component  from 
f«0.10  HZ  to  f«0.05  HZ  and  larger  wave  angle  within  60 
seconds  after  the  monochromatic  wave  has  reached  the 
steady  state.  As  shown  in  Figure  4.15,  the  evolution 
of  the  energy  spectrum  in  two-dimension  is  similar  to 
that  in  one-dimension.  The  variation  of  the  wave 
direction,  moving  away  from  the  center  line  instead  of 
approaching  it,  is  shown  in  Figure  4.16. 


CHAPTER  5 


RANDOM  NAVE  MODEL  FOR 

NON-STATIONARY  WAVE  SPECTRAL  TRANSFORMATION 


In  this  chapter,  discussion  is  focused  on  the 
random  wave  model  which  consists  of  three  sub-models  : 
1)  swell  wave  model,  2)  wind  wtave  model,  3)  hybrid  wave 
model.  The  field  data  comparison  is  also  included.  A 
schematic  view  of  ocean  environment  which  Is  the  same 
as  Figure  3.3  is  illustrated  in  Figure  S.I.  To  apply 
the  numerical  scheme,  A x*5000m,  Ay*5000ra  and  At*60 
seconds  are  used  in  the  three  sub-models. 


5.1  Swell  Wave  Model 


The  swell  wave  model,  without  the  consideration 
of  local  wind  generation,  computes  wave  transfromation 
when  the  offshore  boundary  condition  is  unsteady.  In 
the  example  given  here,  Pierson  -  Moskowitz  (P-M) 
spectrum,  a  fully  developed  sea  of  the  form. 
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1j  chos«n  as  tha  Initial  condition.  The  offshore 
boundary  is  treated  as  open  boundary  vhlch  is  affected 
by  the  offshore  wind  variation. 

Two  different  offshore  open  boundary  conditions 
are  illustrated:  one  for  increasing  wind  speed,  hence, 

growing  spectrum  and  the  other  for  decreasing  wind. 

% 

First,  the  decreasing  wind  case  is  examined.  The 
initial  condition  is  established  by  inputting  a  wave 
spectrum  at  offshore  boundary  (location  6II).  This 
spectrum  corresponds  to  a  fully  developed  sea  of  7 
m/sac  wind  and  a  main  direction  20  degree  normal  to  the 
shore.  The  steady  state  spectra  at  boundary  (location 
6II),  location  SII  and  location  2II  are  denoted  as  t*0 
min  curves  in  Figure  S.2.  The  wave  condition  at  the 
offshore  is  then  allowed  to  vary  by  diminishing  the 
wind  speed  from  7  m/sec  to  4  m/sec  in  one  hour  at 
constant  rate  with  the  spectral  curves  marked 
successively  as  20,40,60  (min)  from  time  0.  The 
corresponding  wave  conditions  at  location  SII  and 
location  2II  for  the  successive  times  are  as 


116 


illustrated.  The  affect  of  non-s tat ionary  vave 
transformation  is  clearly  revealed  here.  By  comparing 
spectra  at  boundary  (location  6II)  and  at  location  2II, 
it  can  be  seen  that  although  the  wave  conditions 
changed  drastically  at  offshore  within  one  hour,  the 
waves  in  nearshore  region  were  hardly  affected  during 
this  time  period. 

Figure  5.3  provides  another  example  with  a 
situation  just  opposite  to  the  previous  one  in  that  the 
initial  wind  condition  is  4  m/sec  and  begins  to  pick  up 
in  an  hour  at  steady  rate  to  a  final  speed  of  7  m/sec. 

Field  data  measured  by  Ijima  (1953)  is  shown  in 
Figure  5.4  which  illustrates  the  time  change  of  wave 
spectrum  during  a  decreasing  period  of 
typhoon-generaced  swell.  The  effect  that  the  wave 
components  disappear  as  wind  speed  decreases  is 
obs  erved . 

5.2  Wind  Wave  Model 

The  wind  wave  model  is  concerned  with  the 
effect  of  local  wind  generation.  A.  calm  sea  is  chosen 
as  the  initial  condition  with  an  open  offshore 
boundary. 
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Time  Changes  of  Wave  Energy  Spectrum  in  the 
Decrease  of  Swell  of  a  Typhoon  (After  Ijima,  1957) 
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Th«  vind  beginning  to  build  up  from  0  to  7 
m/sec  in  the  first  half  hour  and  then  maintaining  7 
m/sec  for  the  second  half  hour,  blows  over  a  calm  sea 
in  the  direction  perpendicular  to  the  shoreline.  As 
shown  in  Figure  S.S,  the  spectral  peak  shifts  tovrard 
lower  frequency  with  increasing  time.  This  phenomenon 
is  consistent  with  the  field  observation  in  the  ocean. 

5.3  Hybrid  Wave  Model 

Based  on  linear  superposition,  the  swell  wave 
model  and  the  wind  wave  model  can  be  combined  into  a 
hybrid  wave  model  to  predict  the  real  ocean  waves . 

Here,  the  offshore  wave  condition  is  the  same 
as  given  in  the  first  case  of  the  swell  wave  model, 
i.e., waves  at  outer  boundary  diminish  from  a  spectrum 
corresponding  to  7  m/sec  wind  to  that  of  a  4  m/sec  wind 
in  one  hour.  In  the  region  of  consideration,  a  local 
wind  of  the  same  nature  as  given  in  the  wind  wave 
model,  i.e.,  wind  begins  to  build  up  from  0  to  7  m/sec 
in  the  first  half  hour  and  then  maintains  at  7  m/sec 
for  the  second  half  hour,  is  also  occurring 
concurrently.  The  corresponding  wave  conditions  at 
location  5II  and  location  2II  are  shown  in  Figure  5.6. 
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iriation  o£  Energy  Spectrum  due  to  Decrease  of 
^ell  and  Increase  of  Wind  V^ave  at  Location  5II 
Location  2X1 
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The  ipcctra  are  now  dual-peaked;  the  peak  due  to  swell 
shifts  toward  higher  frequencies  as  the  offshore  wave 
diminishes  in  strength  whereas  the  peak  corresponding 
to  local  wind  wave  shifts  towards  the  lower  frequencies 
as  the  local  wind  picks  up.  Eventually,  the  wave 
spectrum  of  the  local  wind  wave  will  dominate  the  whole 
wave  spectrum.  This  phenomenon  is  very  clearly  shown 
in  Figure  S.7  from  Long  and  Hasselmann  (1979). 

5.4  Comparison  with  the  North  Sea  Field  Data 

Documented  field  data  with  sufficient 
information  for  the  purpose  of  comparing  with  the 
present  numerical  model  is  scarce.  The  recent  MARSEN  I 
experiment  was  conducted  during  the  months  of  September 
and  October  1979.  Waves  were  measured  in  shallow  water 
by  a  series  of  wave  gages  as  well  as  in  deep  water  by  a 
number  of  pitch-and-roll  buoys .  It  is  an  attempt  to 
understand  the  characteristics  of  wave  spectral 
transformation  from  deep  to  shallow  water. 

The  experiments  were  performed  in  the  German 
Sight  around  the  Island  of  Helgoland  and  Sylt.  Figure 
5.3  shows  the  locations  of  the  nearshore  site  and 
deepwater  measurements  in  the  MARSEN  test  region. 


Figure  5 . 
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Frequency  Spectrum  Indicating  Wind  Sea  and  Swell 
Peaks  (After  Z<ong  and  Hasselmann,  1979) 


Moya  and  ataeiona  In  tha  Caraan  Dl9ht 
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Deepwater  wave  information  wa3  obtained  from 
the  pitch-and-roll  buoy  maintained  by  Deutsches 
Hydrographis ches  Institute  of  Germany.  The  shallow 
water  measurements  by  Wang  et  al.  (1980)  with  emphasis 
on  surf  zone  characteristics  were  carried  out  at  the 
Island  of  Sylt.  The  distance  between  the 
pitch-and-roll  buoy  and  the  nearshore  site  is 
approximately  70  km. 

Figure  5.9  shows  the  instrument  arrangement  in 
the  nearshore  zone.  The  wave  gauges  at  locations  1225 
m  (approximately  1100  m  from  shoreline  ,  1225  m  denotes 
the  distance  in  meter  from  a  base  line  which  is 
approximate  100  m  inland  frcxn  shoreline) ,  940  m,  and 
225  m  are  bottom  mounted  echo-transceivers  as  developed 
by  Fuhrentholz  Laboratory.  These  echo  sounders  are 
capable  of  measureing  water  surface  variation  to  ±  0.5 
cm  for  water  depths  up  to  90  m.  The  shallow  water 
gauges  are  staff-mounted  pressure  transducers  of  type 
MDS  76  as  manufactured  by  H.  Maihak  AG,  Hamburg  with  a 
pressure  range  of  0  -  1  kg/ cm*.  The  current  meters  are 
COMEX  e.’ ectromagnetic  two  component  type. 

Data  were  taken  with  recording  length  of 
approximately  20  minutes.  Spectral  analysis  was 


performed  using  the  FFT  technique  developed  by  Cooley 
and  Tukey  (1965)  with  Hanning  Window  (Blackman  and 


Tukey,  1958)  for  data  smoothing.  All  the  data  sets 
were  sampled  uniformly  at  a  time  interval  i^t«0.5  sec  to 
a  total  data  point  N«2048  or  1024.  In  order  to 
maintain  good  resolution,  a  degree  of  freedom  DF«40  for 
N«2048  and  OFa20  for  H«1024  was  selected  as  the 
optimum. 

The  mean  water  depth  at  location  of 
pitch'-and-roll  buoy  and  at  location  1225  ro  are  23  m  and 
8  m,  respectively.  The  wave  data  were  still  being 
analyzed  and  sorted  out.  So  far,  we  have  indentified 
only  one  set  of  data  that  marginally  suitable  for  the 
present  purpose  (simultaneous  wave  recordings  with  all 
the  necessary  input  information  documented) .  The 
results  of  this  data  set  which  was  taken  on  Sept.  24, 
1979  are  reported  here.  First  of  all,  the  wind 
condition  of  this  particular  day  is  shown  in  Figure 
5.10.  The  wind  speed  variation  was  rather  moderate  for 
the  whole  day  a<maewhere  around  3  -10  m/sec  for  the  bulk 
part  of  the  day;  the  wind  direction  was  from  SWS 
during  0  to  3  o'clock  local  time  and  then  suddenly 
shifted  to  H  from  3  o'clock  till  midnight.  Since  the 


wind  direction  changes  drastically  at  9  o'clock,  wind 
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wave  generated  from  that  time  meets  the  criteria  of  the 
local  wind  model  described  above. 

Wave  data  were  collected  at  offshore  location 
from  8  o'clock  at  approximatly  2  hour  interval.  The 
energy  spectra  and  their  corresponding  directional 
information  are  given  in  Figure  5.11.  Shallow  water 
waves  were  recorded  at  location  1225  m,  during  the  same 
period.  A  rough  calculation  indicated  that  the 
dominant  wave  component  (  0.16HZ)  would  take 
approximately  2  hours  to  travel  from  the  offshore  site 
to  the  nearshore  zone.  Thus,  to  predict  wave  condition 
at  the  nearshore  site,  the  input  at  offshore  location 
should  be  chosen  about  2  hours  in  advance.  Based  upon 
this  approximation,  deepwater  waves  at  12:14  pm  were 
used  as  input  to  compute  ths  condition  at  location  1225 
m  at  14:00  pm.  The  bottom  friction  coefficient  is 
chosen  as  0.01  because  of  the  sandy  bottom.  The  grid 
system,  with  Az«5000  m,  <&y»5000  m  and  i^t^SO  seconds 
which  used  in  the  computation  is  shown  in  Figure  5.12. 

The  comparison  of  numerical  results  and  field 
data  is  shown  in  Figure  5.13  and  Figure  5.14.  The 
contributions  to  the  total  energy  spectrum  from  local 
wind  generation  and  swell  transformation  are  indicated 
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Figure  5.12  Grid  System  used  in  the  Computation 


133 


2  2  2  2 

o  o  o  o  o  • 

(sas-SRBMMiZSNsa  nwiaads 


Figure  5.14  Comparison  between  Hybrid  Wave  Model  and  Field 
Data  at  14t00  at  Location  1225ro 
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in  Figure  5.13.  As  can  be  seen,  for  the  specific  case, 
local  wind  waves  were  the  major  contributor.  The  swell 
component  while  important  to  the  peak  energy  was  minor 
as  a  whole.  A  majority  of  the  wave  components  present 
in  the  offshore  region  actually  never  arrived  at  the 
nearshore  zone.  As  we  have  stated  earlier,  this  set  of 
data  is  really  marginal  for  verification  purpose 
because  of  this  minor  contribution  due  to  swell 
components . 


..1  111.  4  ..  I  HJ 


CHAPTER  6 


SUMMARY  AND  CONCLUSIONS 


A  mathematical  model  which  can  predict  the 
stationary  and  non- stationary  wave  spectral 
transformation  in  shallow  water  of  irregular  bottom  has 
been  developed.  The  model  contains  the  physical 
dynamic  processes  which  are  conserved,  such  as  wave 
refraction  and  shoaling,  as  well  as  those  which  are  not 
conserved,  such  as  wave  energy  generation  and 
dissipation.  The  numerical  results  can  be  used  to 
guide  further  research  in  wave  phenomenon  of  which  very 
little  is  known. 

The  theory  of  wave  spectral  transformation  is 
presented  in  Chapter  2.  Instead  of  using  the 
parametrical  method  which  was  developed  by  Hasselmann 
(1976),  the  theory  is  based  on  the  physical  method. 
Both  methods,  physical  and  parametric,  have  their 
relative  advantages  and  limitations.  The  problem 
becomes  much  more  complicated  when  the  shallow  water 
effect  is  to  be  considered.  Several  additional 
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effects,  such  as  wave  shoaling,  refraction  and  bottom 
dissipation  must  be  taken  into  account.  Furthermore, 
the  non-linear  wave-wave  interactions  become 
increasingly  important  in  shallow  water.  The  physical 
method  includes  the  former  effects  but  excludes  the 
latter  effect,  whereas,  the  parametrical  method  can 
only  handle  the  local  wind  fields  which  are 
characterized  by  extremely  large  spatial  and  temporal 
gradients.  Based  on  the  reasons  specified  above,  the 
physical  method  which  can  handle  the  time  varying  wind 
condition  in  shallow  water  is  chosen  as  the  basis  of 
wave  spectral  transformation.  The  results  are  found  to 
be  reasonable. 

The  theory  consists  of  three  kinematic 
conditions,  which  handle  the  wave  refraction  ,  and  one 
dynamic  condition,  which  controls  the  wave  shoaling, 
wave  generation  and  wave  dissipation.  In  the  wave 
generation,  two  mechanisms  are  considered:  one  is 
Phillips'  resonant  mechanism  which  causes  the  initial 
wave  growth  and  the  other  is  Miles'  instability 
mechanism  which  controls  the  major  wave  growth.  In  the 
wave  dissipation,  there  are  also  two  mechanisms.  The 
bottom  friction  and  wave  breaking  are  chosen  from 
Hasselmann  and  Collins'  derivation  and  Kitaigorodskii ' s 
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shallow  water  equilibrium  energy  spectrum, 
respectively . 

The  numerical  scheme  of  stationary  wave 
spectral  transformation  is  analyzed  in  Chapter  3. 
Based  on  Noda ' s  numerical  procedures,  the 
non-homogeneous  terms  which  are  the  effect  of  the  wind 
generation  and  bottom  friction  are  also  added  to  be 
solved.  Shoaling  and  refraction  are  well  handled  in 
the  model.  In  low  fequency  components,  the  influence 
of  the  bottom  friction  is  found  to  be  important  in 
shallow  water  since  the  long  wave  component  has  already 
felt  the  bottom.  However,  in  high  frequency 

components,  energy  generation  or  transfer  are  more 
important  than  the  bottom  friction.  In  the  steady 
local  wind  generation,  the  short  wave  component  exceeds 
the  maximum  wave  steepness  and  follow  the  breaking 
criterion  of  Kitaigorodskii ' s  shallow  water  equilibrium 
energy  spectrum.  Whereas  the  long  wave  component  has 
not  reached  the  maximum  wave  steepness  and  continues  to 
grow  until  the  total  energy  spectrum  reaches  the 
depth-limited  wave  condition. 

The  numerical  scheme  of  non-stationary  wave 
spectral  transformation  is  analyzed  in  chapter  4.  The 
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problem  is  solved  by  using  the  two-step  Lax-Wendroff 
scheme.  Numerical  verification  is  based  on  an  example 
which  has  the  same  type  as  the  energy  transport 
equation.  The  exact  solution  and  the  approximate 
solution  agree  very  well  (no  more  than  2%).  It  shows 
that  the  present  numerical  scheme  can  handle  both 
transient  equations  (one  is  the  wave  conservation 
equation  and  the  other  is  the  energy  transport 
equation) .  The  characteristics  of  the  wave  spectral 
transformation,  in  which  the  variation  of  wave 
component  at  the  boundary  affects  far  locations,  are 
shown  in  the  example  of  monochromatic  wave  in  one  and 
two  dimensions.  The  effect  which  the  propagation 
velocity  of  long  wave  is  faster  than  that  of  short  wave 
is  clearly  revealed. 

Based  on  the  linear  assumption,  the  random  wave 
model  for  non-stationary  wave  spectral  transformation 
is  discussed  in  Chapter  5.  The  problem  is  handled  by 
developing  three  separate  models;  1)  swell  wave  model 
handles  the  offshore  wind  variation  without  local  wind 
generation,  2)  wind  wave  model  handles  the  local  wind 
generation,  3)  hybrid  wave  model  linearly  combines  the 
swell  wave  model  and  wind  wave  model  to  predict  a  real 
ocean  wind  waves.  The  numerical  results  show 
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reasonably  agreement  with  the  North  Sea  field  data 
except  the  wave  components  at  high  and  low  frequency. 
This  is  likely  to  be  a  consequence  of  the  non-linear 
wave-wave  interactions  which  is  not  considered  in  this 
model  and  of  the  assumed  constancy  of  the 
Kitaigorodskii ' s  "  B  "  constant  for  the  equilibrium 
condition. 

Wave  spectral  transformation  in  shallow  water 
is  important  for  the  study  of  analyzing  and  predicting 
wave  conditions  in  nearshore  zone  due  to  hurricanes  or 
storms  sweeping  through  offshore  region.  In  the 
present  model,  only  the  linear  problem  is  considered. 
Ideally,  the  non-linear  wave-wave  interaction  terms 
should  be  included.  But,  it  has  not  been  done  so  for 
two  practical  reasons  as  follows: 

1)  Lack  of  suitable  model.  Phillips  (1960)  made  the 
first  investigation  of  nonlinear  resonant  interactions 
among  particular  groups  of  wave  components. 
Hasselmann,  in  a  series  of  three  papers  (1962,1963a  and 
1963b)  studied  the  rate  of  energy  exchange  in  a  random 
wave  field.  He  assumed  that  the  probability  density 
distribution  of  the  surface  displacement  is  normal  and 
calculated  the  net  rate  of  energy  exchange  from  one 
wave  to  its  wave  components  in  the  resonant  tetrad. 
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Unfortunately,  the  results  he  got  are  neither 
mathematically  smooth  nor  stable.  Barnett  (1968) 
modified  Hasselmann's  parameterized  function  in  a 
manner  that  will  increase  the  calculation  speed  but  the 
result  is  only  an  approximation.  Based  on  the  narrow 
band  spectrum,  Longuet-Higgins  (1976)  and  Fox  (1976) 
applied  the  evolution  of  wave  packets  in  three 
dimensions  to  the  Hasselmann's  resonant  energy  transfer 
result.  The  simplification  improved  calculation  speed 
and  stability.  All  the  above  studies  are  restricted  to 
deepwater  waves.  Herterich  and  Hasselmann  (1960) 
modified  Longuet-Higgins  results  for  the  computation  of 
the  rate  of  wave  energy  transfer  in  finite  depth  water 
for  a  narrow  band  spectnam.  However,  their  results 
cannot  be  applied  to  shallow  water  areas  where  the  kh 
value  is  between  0.3  and  0.7.  Resio  (1981)  modified 
Hasselmann* s( 1976 )  parameterized  function  to  be  a  new 
parameterization.  The  result  is  still  not  applicable 
because  it  is  restricted  to  deepwater  waves. 

2)  Lengthy  computation.  The  numerical  program  would 
have  become  unwieldily  complicated  at  this  stage  that 
would  submerge  us  in  numerical  manipulation. 

Actually,  it  is  a  drawback  to  omit  the 
non-linear  wave-wave  interaction  terms  which  by  many 


rascarchcrs  ar«  regardad  as  an  important  paramatar  in 
shallow  watar.  Tharafora,  a  mora  axtansiva  study 
should  account  for  tha  a f facts  of  non-linaar  wava-wava 
intarations . 
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